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Abstract 


The  dynamic  coUiq>se  behavior  of  thin  composite  cylindrical  shells  subjected  to 
transverse  point  loads  is  considered.  The  dynamic  behavior  of  the  undamped  cylindrical 
shells  of  varying  radii  is  analyzed  with  a  finite  element  formulation  that  incmporates  all 
nonhnear  Green’s  strain  terms  in  tte  in-plane  directions;  transverse  shear  strains  are 
liiwar  and  vary  parabolically  through  the  shell  thickness.  This  formulation  is  denoted  as 
the  simplified  large  displacement/rotation  (SLR)  theory.  Graphical  representations  are 
used  extensively  to  examine  dynamic  behaviors  not  noted  in  similar  research  that  utilized 
the  DSHELL  tinite  element  code.  Composite  plate  behavior  is  also  explored. 
Comparison  of  the  higher-order  formulations  to  Donnell  and  von  Karman  models, 
modified  with  transverse  shear  flexibility,  is  also  made.  Several  deep  shell  configurations, 
using  the  SLR  and  Donnell  models,  are  analyzed  in  an  effort  to  determine  the  maximum 
displacement  and  rotational  limits  of  these  formulations.  The  nonlinear  features  of  the 
plates  and  shells  imply  the  potential  for  chaotic  behavior.  Various  techniques  are  used 
to  characterize  the  chaotic  behavior  of  these  undamped  sheUs  in  pre-  and  post-collapsed 
states. 
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Comparison  of  the  Dynamic  Behavior  of  Composite  Plates  and  Shells  Incoiporating 
Green’s  Strain  Terms  with  the  von  Karman  and  Donnell  Models 

I.  INTRODUCTION 

There  has  been  considerable  effort  directed  towards  the  study  of  composite  materials 
because  of  their  advantages  over  conventional  isotropic  materials.  The  high  strength-to- 
weight  properties  of  many  composites  are  ideally  suited  towards  aerospace  i4>plications. 
Composites  are  typically  formed  in  lamina  layers  with  the  high  strength  properties  usually 
oriented  in  the  in-plane  directions.  Thus,  composites  can  be  made  lighter  and  thinner  than 
conventional  isotropic  materials  in  a  given  ^plication.  However,  composites  may  be 
substantially  weaker  in  the  transverse  directions,  especially  if  the  composite  is  thin. 
Therefore,  depending  on  the  application,  the  stability  of  a  composite  may  be  an  important 
consideration. 

The  technical  literature  is  rich  in  research^’  efforts  to  model  the  behavior  of 
composite  structures.  Since  the  strength  properties  in  a  lamina  may  vary  drastically  from 
ply  to  ply  with  respect  to  a  given  reference,  composites  are  clearly  more  complicated  to 
analyze  than  isotropic  materials.  In  an  attempt  to  reduce  the  complexity  of  the 
mathematical  model,  analysts  typically  make  simplifying  assumptions  especially  if  the 
laminate  is  thin.  The  model  formulations  in  this  research  exploit  this  dimensionality  as 
well.  When  applications  involve  small  deflections  and  small  rotations,  many  classical 
theories  neglect  or  simplify  nonlinear  tenns  and  yet  retain  good  accuracy.  However, 
when  deflections  become  many  tunes  larg^  than  the  laminate  thickness,  a  more 
thoughtful  approach  towards  including  higher-order  nonlinear  terms  is  required. 
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Linear  bifurcation  analysis  is  not  suitable  in  predicting  the  coll{^>se  behavior  of  shells. 
For  small  deflections  and  rotations,  linear  bifurcation  analysis  may  adequately  predict  the 
nonlinear  stability  limits  of  some  structures,  although  some  degree  of  error  is  expected 
(11:466).  This  is  generally  true  for  Euler  columns  and  axially  loaded  flat  plates  since 
these  structures  experience  little  stress  distribution  or  geometry  changes  prior  to  the  onset 
of  the  buckling  instability.  On  the  other  hand,  thin  composite  shells  may  not  collapse 
until  moderate  deflections,  rotations,  and  stress  redistributions  have  occurred.  Therefore, 
it  is  expected  that  a  real  structure  may  collapse  at  a  load  quite  different  than  that 
predicted  by  a  linear  bifurcation  buckling  analysis. 

Von  Karman  is  considered  the  first  person  to  reconcile  the  buckling  of  cylindrical 
shells,  observed  in  experiments,  with  a  theoretical  formulation  which  incorporates  large 
deflection.  "The  main  feature  of  the  large  deflection  q)proach  is  that  one  or  more 
nonlinear  partial  differential  equations  are  solved  approximately.  The  post-coU^se  curve 
branching  from  the  fundamental  equilibrium  state  is  subsequently  detennined"  (11:465). 
In  his  woik  with  cylindrical  shells  subjected  to  axial  loads  and  external  pressure,  von 
Karman  showed  that,  at  a  load  or  pressure  level  less  than  the  value  predicted  by  linear 
bifurcation  analysis,  the  deforming  cylindrical  shell  jumps  (i.e.,  snaps)  into  another 
equilibrium  state,  corresponding  to  the  coU^sed  configuration.  A  similar  {q)proach  was 
taken  by  Dennis  and  Palazotto  (14);  they  also  observed  similar  .ollapse  characteristics 
in  their  work  with  transversely  loaded  shells.  Their  work  is  the  basis  for  the  tindings 
contained  in  this  research. 

In  recent  years,  a  new  scientific  approach  to  analyzing  complicated  and  seemingly 
unpredictable  behavior  has  developed: 
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This  new  paradigm  goes  by  the  name  of  Chaos,  referring  to  a  scientiric 
philosophy,  an  approach,  and  a  set  of  methods  to  deal  with  the  manifestations  of 
chaos  in  the  physical  sciences.  To  a  large  extent  the  enthusiasm  that  has 
developed  for  Chaos  is  the  result  of  the  breadth  of  its  applications.  These 
£^)plications  for  Chaos  for  understanding  complex  and  ui^redictable  behavior 
range  across  the  spectrum  of  scientitic  disciplines.  Indeed,  Chat^  has  much  the 
same  flavor  as  classical  thermodynamics  in  that  the  fundamental  ideas  and  results 
seem  applicable  to  a  wide  variety  of  different  physical  systems.  There  is  probably 
no  physical  system  exhibiting  unpredictable  behavior  that  is  presently  not  being 
scrutinized  through  the  lens  of  Chaos  by  someone  (16:vii) 

By  definition  chaotic  vibrations  arise  from  physical  systems  whose  input  forces  or 
excitations  are  known  (determirristic).  However,  knowledge  of  all  influences  acting  on 
a  system  and  the  exact  laws  governing  the  behavior  does  not  guarantee  accurate 
prediction  of  how  a  system’s  response  will  evolve  over  time.  Even  with  huge  advances 
in  computer  technology,  Peitgen,  Jflrgen,  and  Saupe  (15:viii)  assert  that  the  ability  to  see 
more  accurately  into  future  developments  is  unattainable.  Moreover,  they  maintain  that 
stricter  determinism  and  apparently  accidental  development  are  not  mutually  exclusive, 
but  rather  that  their  coexistence  is  more  the  rule  in  nature.  Therefore,  as  a  researcher 
attempts  to  accurately  model,  to  a  high  degree,  the  behavior  of  a  physical  system,  some 
of  the  phenomenon  that  he  or  she  may  observe  might  be  best  analyzed  with  the  tools 
developed  from  this  new  scientific  approach. 


1.1  Historical  Motivation 

The  subject  of  dynamic  buckling  of  shells  has  drawn  the  attention  of  researchers  over 
many  years  because  of  the  numerous  possible  applications.  Applications  for  the  aerospace 
industry  was  one  of  the  primary  motivations  for  the  interest  in  cylindrical  shell  buckling 
due  to  the  possible  performance  and  cost  advantages  that  could  result  from  a  weight- 
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optimized  structure.  The  vast  majority  of  researchers  were  interested  in  the  behavior  of 
cylindrical  shells  when  loaded  with  axial,  radial  pressure,  and  torque  fences  since  these 
combinations  are  the  basic  loading  configurations  for  cylindrical  shell  structures  (5:413). 

For  example.  Brush  and  Morton  (1)  published  one  of  the  first  papers  on  experiments  with 
the  static  buckling  of  shells  with  axial  forces,  with  an  intended  rqrplication  for  missiles 
and  launch  vehicles.  One  advantage  to  the  loading  modes  described  is  that  closed  form 
solutions  may  be  derived,  typically  through  the  use  of  series  solutions.  Flugge  (5:413- 
473)  provides  examples  of  this  technique.  One  of  the  first  American  p{q)ers  that 
examined  the  dynamic  budding  of  cyhndiical  shells  was  written  by  Yao  (28)  and 
examined  the  buckling  of  a  cylinder  under  flexural  loads  and  an  external  step  pressure 
pulse.  The  formulations  used  were  based  on  Mathieu  and  Donnell  type  equations. 

In  the  last  30  years,  there  have  been  several  approaches  and  solution  methods  in  the 
static  and  dynamic  analysis  of  plates  and  shells.  In  order  to  reduce  the  complexity  of  the 
general  elasticity  equations,  many  researchers  utilized  thin  plate  and  shell  assumptions 
which  include  1)  small  rotations  and  deflections,  2)  normals  to  the  shell  datum  surface 
remaining  normal  before  and  after  deformation,  and  3)  transverse  normal  stresses  that  are 
negligible.  A  common  method  is  to  use  a  two  dimensional  formulation  to  predict  plate 
and  shell  behavior.  Noor  and  Burton  (12)  have  classified  the  various  methods  of 
modelling  a  multilayer  plate  with  two-dimensional  models  ({q>plicable  to  shell  theory  as 
well)  into  4  general  jq)proaches.  They  are  the  1)  method  of  hypotheses,  2)  method  of 
expansion,  3)  asymptotic  integration  technique,  and  4)  iterative  methods  and  methods  of 
successive  corrections.  The  first  approach  is  based  on  Kirchoff-Love  formulations  with 
assumptions  on  the  variations  of  displacements  and  strains  (or  stresses)  in  the  thickness 
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direction.  The  second  i^proach  is  based  on  series  expansions  with  respect  to  the 
thickness  parameter  of  the  plate  or  shell.  Power  series,  Legendre  polynomials,  and 
trigonometric  series  have  been  used.  The  third  approach  involves  tte  formulation  of 
recursive  sets  of  two-dimensional  equations  from  the  three-dimensional  elasticity 
equations.  The  fourth  approach  includes  iterative  approximations  of  the  three-dimensional 
elasticity  equations  and  successive  corrections  to  the  two-dimensional  formulations. 

1.2  Development  of  Current  Approach 

Dennis  and  Palazotto  (14)  developed  a  simplified  large  displacementA'otation  (SLR) 
formulation  where  the  in-plane  strain  equations  incorporated  all  nonlinear  Green’s  strain 
terms  while  the  transverse  strain  equations  retained  linear  terms  only.  The  transverse 
strain  equations  satisfy  boundary  conditions  on  the  upper  and  lower  surface  of  the  shell. 
The  kinematics  equations  are  based  on  a  power  series  expansion  with  respect  to  the 
thickness  direction.  The  SLR  formulation  lead  to  the  development  of  the  SHELL 
computer  code  which  was  used  for  the  static  analysis  of  plates  and  sheUs.  This  is  one  of 
the  computer  codes  used  in  this  research.  Comparisons  to  the  previous  woric  of  Sabir  and 
Lock  (14:2(X))  showed  excellent  agreement  between  the  two  approaches.  Sabir  and  Lock 
used  an  assiuned  strain  function  and  the  finite  element  method  to  predict  the  collapse 
behavior  of  cylindrical  shells.  The  SLR  approach  was  also  expanded  to  include  analysis 
of  shells  composed  of  composite  materials. 

Tsai  and  Palazotto  (26)  incorporated  the  usability  to  perform  dynamics  analysis  into 
the  SLR  approach.  Their  DSHELL  computer  code  is  one  other  program  used  in  this 
work.  A  Newton-Rhapson  iterative  method  and  Beta-m  time-marching  integration  scheme 
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is  used  to  solve  the  differential  equations  of  motion.  Comparisons  of  DSHELL  results 
to  the  work  of  Qough  and  Wilson  (3),  whose  arbitrarily  thin  shells  utilized  an  assumed 
constant  transverse  shear,  showed  excellent  agreement.  Greer  (6)  added  the  IDEAS™ 
post-processing  feature  to  the  program  to  gr^hically  display  the  shell  at  any  point  in  the 
loading  history  (this  feature  was  not  used  in  this  work). 

13  Chaos 

An  essentially  two  dimensional  structural  system,  such  as  the  thin  composite  plates 
and  shells,  is  an  extremely  complex  structure  to  analyze,  especially  when  chaotic  behavior 
is  to  be  examined.  Other  work  in  numerical  simulation  of  chaotic  motion  typically 
reduces  the  complexity  by  reducing  the  dimensionality  and  degrees  of  freedom.  These 
simplified  models  still  allow  for  a  diversity  of  nonlinear  features  that  is  necessary  in  a 
truly  chaotic  system.  Researchers  such  as  Siljak  (21),  Moon  (10),  Nayfeh  and  Mook  (12), 
Shaw  and  Holmes  (20),  and  Kounadis  (8)  have  examined  simple  structural  systems 
(usually  with  one  or  two  degrees  of  freedom)  with  geometric  or  material  nonlinearites. 
Simulation  of  nonhnear  behavior  has  also  been  carried  out  using  piecewise-linear 
properties.  For  example,  Thompson  et  al.  (25:291)  examined  the  dynamics  of  a  tower 
using  linear  damping  and  bilinear  stiffness  beams.  Their  results  demonstrated  complex 
dynamic  behavior  that  contained  chaotic  characteristics.  Liaw  (9)  performed  a  more 
extensive  analysis  by  examining  a  one-degree-of-freedom  model  that  had  similar  stiffness 
and  damping  properties.  He  demonstrated  that  the  system  exhibited  chaotic  tendencies 
through  quantitative  (e.g.,  Lyapunov  exponents)  and  qualitative  (e.g.,  bifurcation 
diagrams)  analyses. 
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1.4  Current  Work 


An  overall  purpose  of  this  research  u  to  further  elaborate  on  the  dynamic  behavi(»r 
of  cylindrical  shells  as  predicted  by  the  DSHELL  computer  code.  Thus  far,  only  a  few 
researchers  have  used  the  DSHELL  computer  code  to  examine  the  complex  collapse 
behavior  of  cylindrical  shells.  With  so  many  possible  shell  configurations,  th^  is  the 
potential  to  do  much  more.  Shell  geometry,  lengths,  thicknesses,  boundary  conditions, 
and  various  composite  ply  orientations  and  materials  are  just  a  few  of  the  parameters  that 
can  be  varied. 

One  significant  hurdle  was  recently  overcome  by  Greer's  wodc  (6)  in  defining  a 
method  for  selecting  computational  parameters  that  would  lead  to  nummically  stable 
results.  Early  efforts  with  the  DSHELL  code  lead  to  results  that  were  often 
misinterpreted.  For  example,  it  was  originally  thought  that  DSHELL  could  not  perform 
accurate  post-colli^  analysis  due  to  highly  erratic  results.  However,  Greer  demonstrated 
that  the  prc^r  choice  of  certain  computational  parameters  would  allow  prediction  of 
behavior  long  after  the  shell  had  coUi^sed.  This  ability  to  produce  numerically  stable 
results  facilitated  the  study  of  chaotic  motion  in  the  collapsed  shell. 

In  previous  work,  most  of  the  studies  in  the  dynamic  behavior  of  cylindrical  shells 
centered  on  the  behavior  of  shells  with  shallow  to  moderate  dq^ths.  It  was  not  clear  that 
DSHELL  could  predict  the  behavior  of  deeper  shell  geometries,  nor  was  it  clear  Gieer’c 
method  could  be  q)plied  to  these  geometries,  another  objective  is  to  gain  insight  on  the 
limits  of  the  SLR  models.  Moreover,  there  are  extremely  limited  studies  that  examine 
composite  cylindrical  shell  responses  when  subjected  to  a  sudden  transverse  point  load. 
A  basis  of  comparison  is  therefore  generated  by  using  a  more  linearized  model,  namely. 
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the  modified  £)onnell  formulation.  Finally,  an  attempt  is  made  to  gain  insight  into  the 
chaotic  behavior  of  complex  structures,  such  as  composite  shells.  Towards  these  goals, 
SPARC2,  SPARCIO,  and  Silicon  Gnq>hics  UNIX  computer  systems  were  used 
extensively. 
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//  THEORETICAL  DEVELOPMENT 


2.1  Simplified  Large  DefiectwfU Rotation  (SLR) 

The  finite  element  formulation  used  in  this  research  is  based  on  a  total  Lagrangian 
approach  and  features: 

(1)  all  non-linear  in  plane  Green’s  strain  terms  and  linear  transverse  strain  terms 

(2)  the  capability  for  shell  normals  to  the  midplarre  to  lose  their  mthogmiality  after 
loading. 

(3)  zero  033  normal  stress  due  to  thin  shell  assumptions. 

(4)  transverse  shear  strains  that  vary  parabolically  through  the  shell  thickness  and 
vanish  on  the  top  and  bottom  surfaces  of  die  shell. 

Although  a  shell  is  three  dimensional,  a  two  dimensional  fmmulation  is  used  to  model 
a  thin  shell.  Figure  1  shows  a  typical  shell  element. 


Figure  2.1.  36  DOF  isoparametric  shell  element 


where 


•  ^  =  x/a  and  T|  =  s/b,  the  curvilinear  coordinates  of  the  surface 
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•  u  is  the  displacement  in  the  x  coordinate  direction 

•  V  is  the  displacement  in  the  s  comdinate  directicm 

•  w  is  the  transverse  displacement  in  the  z  coordinate  direction 

•  w,,  is  the  slope  bw/bx  at  the  node 

•  w,^  is  the  sl(^  bwibs  at  the  node 

•  'Fi  is  the  rotation  of  the  normal  to  the  shell  midsurface  in  the  x  directimi  (about 
the  s  axis 

is  the  rotation  of  the  normal  to  die  shell  midsuiface  in  the  j  direction  (about 
thex  axis 

The  general  shell  kinematics  in  orthogonal  curvilinear  coordinates  are 


or 


Ki(5,ti,0=« 

«2(5.’1.C)=v 


l-i.' 

-  ct,  -  c^* 

1  «lj 

1-i’ 

1 

«3(5,11)=  w 


(1) 


«,(5.n,Ol 

«2«.»i.0[= 


«(5.n) 


w(5,ti),,  . 


♦iCS.’i) 


where  k  =  -4/(3h^.  ^  and  Ti  are  in  the  natural  coordinates  and  C  is  the  coordinate  along 
the  mid-smface  normal.  [/?]  is  a  matrix  that  is  a  fimction  of  ^  only.  For  a  cylinder, 
C//?i^.  For  a  flat  plate,  and  ^  0.  /?,  and  R2  are  shell  radii  in  the  C 
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and  coordinate  directions  respectively. 
The  strains  are  defined  as 


(2) 


where  y,j  is  the  Green’s  strain  tensor  and  h,  are  scale  factors  (19:594).  Insertion  of 
equation  (1)  into  equation  (2)  results  in  strain  equations  with  many  linear  and  nonlinear 
trams.  In  this  investigation  only  the  linear  terms  of  strains  e^^  and  £33,  are  retained. 
Palazotto  and  Dennis  (14)  included  transverse  shear  flexibility  in  the  DonneU  and  Von 
Kaiman  strain  equations.  The  in-plane  modified  Donnell  (14:63)  strain  equatirais  for 
cylindrical  shells  are 

‘u  •  »U  *  |'*’j 

'a  =  ♦  I*! 

'66  =  (3) 

and  utilize  the  kinematic  trams  of  equation  (1).  The  modified  Vrai  Karman  plate 
equations  are  similar  except  is  zero.  The  transverse  shear  equations,  after 
linearizations,  are 

«44  =  ^,2  +  tz  +  C*3*(w^  +  tj) 

«S5  =  *  C*3*(w,  +  t,)  (4) 

Shape  functions  are  used  to  relate  the  u(^,Tl)  and  v(^,ti)  to  the  values  of  u  and  v  at 
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each  node  of  the  element. 


v(«.Ti) 

,  M5,»i)4 
♦2(5.^!) 


(5) 


where  [N]  is  a  matrix  of  shi^  fiinctimis.  =  (m  v  w  w.,  w,2  V;  fw >*  1,  2, 

3,  4  and  «  (u  v]j  for >s5,  6,  7,  8.  Quadratic  shi^  fimctitMis  are  used  to  relate 

the  u(^.T\)  and  v(^>il)  to  the  values  of  u  and  v  at  each  node.  Linear  Lagrangian  shape 
functions  are  used  to  relate  and  ^  ^  values  of  u  and  v  at  each  node. 

Hermitian  shape  functions  are  used  to  relate  the  w(^,Tl)  to  the  values  of  w,  w  „  and 
at  each  comer  node.  See  (14:83). 

The  equation  of  motion  is  derived  using  the  HamilUm  princqMd  where  the  variation  in 
the  time  integral  of  the  total  energy  is  set  to  zero: 

a  f\E-T-W;i=  0  (6) 


where  £  is  the  internal  strain  energy,  the  work  of  nonconservative  fences,  and  the  woik 
of  body  forces.  T  is  the  kinetic  energy  and  11^,  is  the  external  woik.  The  virtual  wenk 
for  a  laminated  plate  and  shell  are 
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(8) 


6r=  f  T  f'*{p«v*avj*}</WQ  (9) 

where  iV  =1-3.  Ck-i  ^  Ck  ^  positions  at  the  top  and  bottom  of  the  itth  layer,  fi  is  the 
domain  of  the  mid  surface,  ay®  6.y®,  p‘''\and  c^'^^are  the  stress  tensca, 

virtual  strain  tensor,  body  force  vector,  virtual  displacement  vector,  velocity  vectm*,  mass 
density,  and  damping  coefficient  for  the  ikth  laym*  respectively.  F,  is  the  extmial  force 
vector.  L  denotes  the  number  of  layers  in  the  laminate. 

The  dynamic  equation  of  motion  can  be  derived  from  (6). 


and  are  acceleration,  velocity,  and  displacement  vectors  respectively  at  each 
nodal  point.  P{t)  is  the  external  load  ^plied  to  each  point.  K  can  be  stated  as: 


where  Kq  is  a  constant  stiffness  matrix,  is  a  stiffness  matrix  related  to  linear 
displacement,  and  Aj  is  a  stiffness  matrix  related  to  quadratic  displacement.  The 
consistent  mass  matrix  is  obtained  by  substituting  equation  (1)  into  equation  (9)  and 
integrating  by  parts  to  obtain 

■  LE  /,'■  (i2) 


The  damping  matrix  is  obtained  in  a  similar  manner,  however,  for  this  research,  no 


damping  is  ccoisidered.  However  for  completeness,  damping  terms  will  be  carried 
through  the  develq>ment  below. 

2.2  Numerical  Technique 

Solution  of  the  dynamic  equation  is  accomplished  through  the  Newton-Raphson  and 
Beta>m  methods.  The  Beta-m  method  (7)  is  a  generalization  of  the  Newmark  time 
marching  integration  scheme.  The  Beta-m  method  transforms  equation  (10)  into  a  set  of 
nonlinear  algebraic  equations  where  the  loading  at  each  point  at  a  given  time  is  known. 
A  Newton-Raphson  iteration  method  is  then  used  to  solve  these  equations. 

The  Beta-m  method  has  some  important  advantages  for  numerical  computation.  First, 
it  is  a  single  step  method  and  therefore  requires  only  the  latest  data  from  the  previous 
time  step.  Second,  the  method  is  general;  the  differential  equation  determines  the 
method  order  m  and  the  user  selects  m  integration  parameters  Pq,  P,,  ...,  p^,.  These  Pj 
selections  influence  stability  and  accuracy  of  the  particular  algorithm  chosen  to  solve  the 
differential  equation.  Third,  the  finite  difference  scheme  is  not  required.  The  beta-m 
method  is  defined  by 

=  9*  +  (13) 

where 

«.  - 1  uTh'-^m-k)\,  (w) 

and 

h  =  (15) 

where  A:  =  0,  l,...,m  and  A  is  the  time  increment.  For  this  research,  P^  =  0.5,  P,  =  0.5, 
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and  ^2  1-0- 


Substitution  of  (13)  into  (10)  at  time  results  in 
[b^  +  b,C  *boKiq^  +  V  £/<->)]  A  t;<->  =P.,,  -{Mq^  *Kiq^  +  b^b,  U^’^>)q^  (16) 

Equation  (16)  is  a  set  of  nonlinear  algebraic  equations  which  may  be  solved  by  various 
algorithms.  A  Newton-Rsf>hson  iterative  method  is  used  in  the  DSHELL  code.  At  each 
time  step 

bUiT!  =  bVf^  +  bbUf*  (17) 

where  i  is  the  iteration  number.  Insertion  of  equation  (17)  into  (16)  results  in 

[b^  4  b,C  *b^^q^ + V 

(18) 

where  is  a  tangential  stiffness  matrix  at  deformation  state  which  is 

Kj(U^^  =  Aq  +  (19) 

Equation  (18)  is  solved  by  the  following  algorithm: 

(1)  Given  . at  time  we  seek  solutions  at  r^,. 

(2)  Calculate  qo,  q„...,  from  equation  (14). 

(3)  Given  and  f/n +°i’  fro™  iteration,  we  obtain  the  right  hand  side 
of  equation  (18)  and  update  the  tangential  stiffness  matrix  from  (19). 

(4)  Solve  for  the  unknown,  from  equation  (18). 

(5)  Calculate  the  update  of  solution  vector  from  equation  (17). 
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(6)  Update  . the  (/  +  l)th  iteration  from  equaticxi  (13). 

(7)  If 


for  a  specified  convergence  criteria  e,  then  stq)  (1)  is  repeated  for  the  next  time  step, 
otherwise  step  (3)  is  repeated  until  equation  (20)  is  satisfied.  The  value  of  /  is  the  degree 
of  freedom  (dof)  number  and  L  is  the  total  number  of  dof. 

2.3  Chaos 

"Chaotic  is  a  term  assigned  to  that  class  of  motions  in  deterministic  physical  and 
mathematical  systems  whose  time  history  has  a  sensitive  dependence  on  initial  conditions" 
(10:4).  Moreover,  in  contrast  to  random  motion,  there  are  no  unpredictable  inputs  or 
parameters  in  a  system  that  is  classified  as  chaotic.  Linear  systems  cannot  exhibit  chaotic 
vibrations  (10:39);  the  system  must  have  nonlinear  features.  In  mechanical  systems,  tl^re 
are  several  possible  ways  for  nonlineaiities  to  exist  and  thus  tte  potential  for  chaotic 
behavior  to  occur.  For  example,  nonlinearities  may  exist  in  elastic  elements,  damping, 
material  properties,  geometric  relations,  boundary  conditions,  backlash,  play,  etc.  In  this 
research,  the  mechanical  system  is  elastic  with  nonlinear  geometric  effects  arising  from 
large  displacements  and  moderately  large  rotations. 

Th^  are  several  methodologies  available  to  indicate  the  presence  of  chaotic 
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behavior.  For  this  work,  the  following  methods  are  ctxisidered; 

•  Observation  of  phase-plane  histories 

•  Fourier  spectrum  analysis 

•  Poincar6  mq)s 

•  Calculation  of  Lyapunov  exponents 

The  motion  of  some  given  feature  of  a  mechanical  system  subjected  to  excitation  may 
be  displayed  in  terms  of  a  dimensional  phase-plane  history.  For  example,  a  plot  of 
velocity  versus  position  (a  classic  phase-plane  plot)  of  some  point  on  a 
mechanical  system  rq>resents  a  two-dimensional  map.  If  this  point  on  the  system  were 
in  an  equilibrium  state,  the  phase  plane  history  would  appear  as  a  point;  for  periodic 
motion,  the  trajectory  mt^)  would  appear  as  a  i  losed  curve;  for  quasiperiodic  motion,  the 
trajectory  map  may  be  plotted  as  a  surface  in  three-dimensional  space.  For  chaotic 
systems,  the  phase-plane  history  tends  to  consist  of  orbits  whose  trajectories  tend  to  fill 
up  a  portion  of  the  phase  space.  A  pseudo-phase-plane  mq)  also  yields  similar  results 
for  all  of  the  types  of  systems  discussed.  A  pseudo-phase-plane  map  (also  known  as  the 
embedding  space  method)  plots  q(t  +  At)  versus  q{t)  where  q  is  some  system  variable. 
This  method  is  particularly  useful  when  only  one  variable  is  measured.  The  idea  is  that 
the  signal  q(t  +  At)  is  related  to  its  time  derivative  and  should  exhibit  properties  similar 
to  the  classic  phase  plane  display. 

Another  clue  to  detecting  chaotic  vibrations  is  the  presence  of  a  broad  spectrum  of 
frequencies  in  a  signal  analysis.  The  fast  Fourier  transform  is  a  popular  method  to 
decompose  a  signal  into  a  set  of  sinusoidal  or  harmonic  signals.  It  is  assumed  that  any 
signal  may  be  rq>resented  as  a  function 
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(21) 


where  =  cos  co/  +  i  sin  cor  and  F  is  a  path  along  a  the  complex  plane  co.  However, 
calculation  of  equation  (21)  is  not  efficient.  Fast  Fourim*  transform  algorithms  are 
typically  used  when  data  sets  composed  of  discrete  points  are  to  be  processed.  The  fast 
Fourier  transform  is  defined  by 

m  =  (22) 

i-i 

where  /  and  J  are  integers.  Moon  (10:149)  indicates  that  a  typical  frequency  spectrum 
plot  may  consist  of  many  closely  grouped  frequency  spikes,  over  a  broad  frequency  range, 
that  gradually  diminish  in  strength  as  the  frequency  increases. 

A  quahtative  feature  of  a  phase>plane  mq)  useful  in  the  study  of  chaotic  systems  is 
the  Poincar6  map  cn*  section.  A  Poincar^  map  consists  of  a  sampling  of  points  in  a 
dimensional  map  at  some  predetermined  interval  of  time  or  space.  For  example,  if  one 
establishes  planes  at  even  intervals  of  time  in  a  phase-plane  trajectory  and  m^s  those 
points  (where  the  trajectories  intersect  the  planes)  onto  a  single  plane,  a  Poincard  m^  is 
formed.  For  a  forced  vibration  systems.  Moon  (10:50)  has  recommended  that  for  a  plot 
of  velocity  x^‘^(0  versus  position  x^°’(/„), 

=  nr  +  To  (23) 

where  7  is  the  period  of  the  driving  motion,  Tq  is  a  delay  time  which  is  genm-ally  some 
fraction  of  the  period  T  (See  Wolf  (27)),  and  n  is  an  integer.  For  harmonic  systems,  the 
Poincard  map  would  qjpear  as  distinct  and  finite  points,  the  number  of  points  dependent 
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on  the  parameters  chosen  in  equation  (21).  For  chaotic  systems,  the  Poincar6  map  may 
consist  of  an  arrangement  points  whose  appearance  d^>ends  on  the  type  of  system 
involved.  For  systems  with  little  or  no  damping,  the  Poincar6  map  often  appears  as  a 
cloud  of  unorganized  points  in  the  phase  plane;  Figure  2.2a  is  a  an  example.  For  chaotic 
systems  with  dissipation  involved,  the  map  often  qjpears  as  a  structure  with  "an  infinite 
set  of  highly  organized  points  arranged  in  what  appears  to  be  parallel  lines"  (10:53). 
Figure  2.2b  is  an  example  of  a  Poincar^  map  exhibiting  fractal  patterns  which  are 
associated  with  strange  attractors.  An  attractor  is  a  set  of  points  or  a  subspace  toward 
which  a  time  history  {q)proaches  after  transients  die  out.  For  example,  for  equilibrium  in 
which  no  motion  is  present,  a  Poincar6  mq)  would  be  a  single  point  since  the  phase-plane 
trajectory  would  pierce  all  phase  planes,  in  the  same  location  on  that  plane,  fOT  all  times. 

One  quantitative  measure  of  chaotic  motion  is  the  calculation  of  the  Lyapunov 
exponent.  In  general,  the  initial  conditions  of  a  dynamical  system  may  be  thought  to 
reside  in  an  infinitesimal  sphere  of  diameter  in  a  phase  space.  Chaotic  trajectories  are 
considered  to  map  the  sphere  into  an  ellipsoid  whose  major  axis  changes  according  to: 

d  .  (24) 

where  "k  is  the  Lyapunov  exponent.  For  periodic  or  harmonic  motion,  A.  will  be  less  than 
or  equal  to  0.  For  chaotic  motion,  X  is  positive.  However,  small  amounts  of  data  are  not 
sufficient  to  determine  if  chaotic  motion  is  present.  One  method  averages  several  data 
values  over  a  different  regions  of  the  phase  space.  Wolf  (27:295)  proposed  the  following 
algorithm: 
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(b) 


Figure  2.2  (a)  Poincar6  map  of  chaotic  motion  of  a  buckled  beam  with  low  damping, 
(b)  Poincar6  map  of  chaotic  motion  of  a  buckled  beam  for  higher  damping  showing 
fractal-like  structure  from  a  strange  attractor  (from  Moon  (10:51)). 
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where  N  ccnresponds  to  the  number  of  rq>lacement  stq>s  (into^als)  in  a  set  of  discrete 
data  points.  An  excellent  explanation  of  die  process  along  with  a  ccxnputer  algcnithm  is 
provided  in  the  aforementioned  reference. 
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Ill  PROCEDURES 


In  this  research,  shells  of  various  radii  were  subjected  to  st^wise  transvme  point 
loads.  Figure  3.1  shows  an  example  of  a  cylindrical  shell  and  its  hinged  boundary;  h  is 
the  depth  of  the  shell  and  c  is  the  chord  length  The  curved  edges  (arc  lengths)  are  free 
and  have  a  fixed  length  of  12  inches,  while  the  straight  edges  are  hinged  and  have  a 
length  of  18  inches.  However,  with  the  arc  lengths  fixed  and  radii  prescribed,  the 
curvilinear  shape  of  any  cylindrical  shell  may  be  described  as  a  ratio  of  the  dqith  to  the 
chord  length  as  long  as  the  thicknesss  and  lengths  remain  proportional  to  the  other 
dimensions. 

The  transverse  point  load  is  applied  in  a  stq>-wise  fashion  to  shells  of  various  radii. 
Figure  3.2  illustrates  the  loading  curves  applied  to  the  center  node.  The  initial  loading 
rate  was  held  constant  at  1.4  x  10^  Ibs/s  for  all  cylindrical  shells.  This  loading  rate  was 
the  same  value  used  by  Taylor  (24)  and  Greer  (6).  Cylindrical  stell  radii  of  4,  6,  8,  12, 
20, 60  inches  were  considered.  Figure  3.3  shows  the  various  cylindrical  shell  geometries 
considered.  A  flat  plate  was  simulated  by  increasing  the  radius  to  an  arbitrarily  large 
value  of  900,000  inches.  Figure  3.4  shows  the  flat  plate  mesh.  In  all  cases,  the 
longitudinal  lengths  of  the  composite  is  18  inches.  The  maximum  load  level  for  each 
shell  was  chosen  as  the  static  coUt^se  load  for  that  sheU.  The  SHELL  (14)  computo-  code 
was  used  to  calculate  this  load.  Chien  and  Palazotto  (2)  showed  that  the  non-linear  static 
equilibrium  load  that  causes  buckling  predicts  the  upper  bound  solution  of  the  dynamic 
buckling  load.  The  collapse  loads  were  detennined  by  a  finite  element  method  using  a 
displacement-increment  method  (14:131);  that  is,  displacement  is  specified  and  the 
equilibrium  load  corresponding  to  that  displacement  is  determined  using  Newton- 
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Load, lbs 


Figure  3J.  Loading  curves  applied  to  cylindrical  shells  of  various  radii. 


Radius  =  4.0  in. 


Radius  =  6.0  in. 


Radius  =  20.0  in.  Radius  =  BO.O  in. 


Figure  Meshes  representing  cylindrical  shells. 


Raphson  iterati(Hi.  All  in-plane  nonlinear  terms  and  linear  transvme  shear  terms  are  used 


Figure  3.4.  Mesh  rq>i:esenting  a  flat  plate. 

in  the  formulation  similar  to  the  above.  The  SHELL  [  ]  program  was  used  to  detmmine 
the  static  collapse  loads. 

Due  to  the  extensive  computational  resources  (i.e.,  time  and  memory)  required  to  run 
the  finite  element  code,  all  shell  configurations  were  limited  to  a  7x7  quarter  element 
Taylor  (24)  mesh.  When  compared  to  finer  Silva  (22)  meshes,  the  Taylor  meshes  fell 
within  5%  of  the  criteria  set  forth  by  Silva  based  on  peak  loadings;  thus  the  coarser  mesh 
was  deemed  adequate  to  ct^ture  the  gemral  behavior  of  the  shell  as  it  colliq>ses. 
Moreover,  of  the  many  possible  degrees-of-freedom  that  could  be  analyzed,  it  is  assumed 
that  the  generalized  displacement  associated  with  the  degree-of-fieedom  of  the  load 
dominates  the  behavior  of  the  system.  Therefore,  the  center  node  displacements  are  of 
primary  interest  throughout  this  woik. 

The  shell  prq)erties  were  modeled  after  a  24'ply,  0. 12  inch  thick,  Hercules  AS4-3S01-6 
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[(V906],  graphite  epoxy.  Hie  material  properties  are  listed  in  Table  3.1.  The  shell  is 
assumed  to  behave  in  a  perfectly  elastic  fashicxi  with  no  interlaminate  failure  and  no 
damping.  The  density  is  0.00015088  slugs^^. 

In  a  dynamically  loaded  shell,  the  choice  of  time  step  and  ccmvergence  tolerance 
is  a  crucial  parameter  for  producing  numerically  stable  results  fen*  the  integrating  scheme 
used  in  this  research.  In  particular,  the  time  step  required  to  c^ture  the  behavior  of  the 
shell  after  budding  must  be  significantly  smaller  than  the  time  step  before  collapse.  For 
a  12  inch  radius  shell,  Greer  (6)  determined  that  the  time  step  required  to  capture  the 
prebuckling  behavior  of  the  shell  is  i^yproximately  x,/20,  where  t,  is  the  period  (in 
seconds)  of  the  first  fundamental  frequency  determined  from  an  eigenvalue  analysis. 
Using  this  time  step  for  post  colli4)se  analysis  results  in  numerically  unstable  results. 
Numerical  instability  typically  leads  to  erratic  data  or  unbounded  motion.  The  time  stq) 
required  for  numerically  stable  results  for  the  12  in.  radius  shell  was  recommended  to  be 
T^lbO,  where  Tj  is  the  period  of  the  highest  frequency  experienced  by  a  free  edge  node 
as  the  shell  oscillates  with  an  tilled  load  that  does  not  result  in  collapse;  the  value  of 
Tj  was  determined  by  a  fast  Fourier  transform  analysis  of  the  displacement  data  at  the 
free  edge  node. 


Table  3.1.  Physical  pr(^)eities  of  AS4-3501-6  gr^hite  epoxy  composite  (psi) 


E, 

E, 

Gn 

Gi3 

G23 

HH 

18.844(10®) 

1.468(10®) 

0.91(10®) 

0.91(10^*) 

0.45(10^) 

0.26 

In  Greer’s  work,  the  value  of  Xj  corresponded  closely  to  the  second  natural  frequency 
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period  of  the  12  in.  radius  cylindrical  shell.  Therefore,  in  an  attempt  to  reduce  the 
wQildoad,  the  value  of  Tj  in  this  research  was  initially  chosen  as  the  naural  preiod 
corresponding  to  the  second  eigenvalue.  The  value  of  was  then  rounded  down  to 
the  nearest  single  significant  digit.  The  DSHELL  (26)  computer  code  was  used  to 
determine  the  eigenvalues.  Table  3.2  lists  the  first  three  natural  frequencies. 

Table  3.2.  First  three  natural  frequencies  of  composite  cylindrical  shells  of  varying  radii 


Radius 

(in.) 

1st  Frequency  (Hz) 

2nd  Frequency  (Hz) 

3rd  Frequency  (Hz) 

4 

460 

1339 

1535 

6 

486 

1070 

1524 

8 

497 

911 

1520 

12 

505 

757 

1512 

20 

508 

655 

1251 

60 

479 

547 

699 

TECPLOT™  was  used  to  display  mesh  geometries.  The  output  of  the  DSHELL 
program  produces  changes  in  the  degree  of  freedom  variables  with  respect  to  the  original 
curvilinear  coordinate  system.  Therefore,  the  new  position  may  be  found  by  algebraically 
adding  the  changes.  In  order  for  TECPLOT^  to  grq)h  the  results,  a  transformation  from 
the  curvilinear  coordinate  system  to  a  Cartesian  system  is  first  required.  For  simplicity, 
only  the  displacement  coordinates  for  the  comer  nodes  (mid-side  nodes  have  in-plane 
displacements  only)  were  utilized  for  graphical  display  purposes.  The  transformation 
requires  simple  trigonometric  relations.  Simple  computer  algorithms  were  written  and 
used  to  carry  out  the  transformation. 
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IV  RESULTS  AND  DISCUSSION 


There  are  several  objectives  addressed  in  this  research.  A  major  objective  is  to 
observe  any  dynamic  cylindrical  shell  behavior  characteristics  not  noted  in  other  similar 
studies  with  the  DSHELL  program.  Anotl^r  objective  is  to  explore  displacement  and 
rotational  limits  under  which  the  SLR  and  the  modified  Donnell  may  be  i^ipUcable  in  a 
dynamics  application  for  plate  and  cylindrical  shell  structures.  Comparisons  between  the 
results  predicted  by  the  SLR  and  Donnell  formulations  are  to  be  made.  Numerical 
stability  characteristics  are  also  examined  as  the  shell  radius  varies.  One  other  goal  is  to 
demonstrate  whether  or  not  these  point-loaded  cylindrical  shells  exhibit  chaotic  motions 
at  certain  phases  in  the  loading  history. 

Previous  research  has  indicated  that  the  SLR  formulation  for  static  analysis  produces 
acceptable  results  when  no  more  than  approximately  15  percent  of  a  shell  structure 
undergoes  up  to  40  degrees  of  rotation  (14).  For  these  shells,  maximum  rotation  would 
be  expected  to  occur  at  the  hinged  boundaries.  With  these  particular  shell  geometries,  the 
8  in.  radius  shell  would  be  expected  to  be  at  or  slightly  above  the  limits  of  the  maximum 
rotation  after  collapse  if  a  static  analysis  was  applied.  Therefore,  of  the  shells  examined 
in  this  research,  the  behavior  characteristics  of  the  8  in.  radius  cylindrical  shell  are 
consider  first. 

4.1  Eight  in.  Radius  Cylindrical  Shell 

Figure  4.1  shows  the  center  node  displacement  of  the  8  in.  radius  cylindrical  shell. 
Upon  initial  ^plication  of  the  load,  both  formulations  show  that  there  is  an  initial  and 
nearly  constant  velocity  increase.  After  t^proximately  1  in.  of  displacement,  the  center 
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Figure  4.1.  Center  node  displacements  of  composite  cylindrical  shells  subject  yd  to  a 
transverse  step  load  tq^plied  to  the  center.  Radius  =  8. 


node  of  the  SLR  shell  reduces  its  average  velocity  as  the  center  node  oscillates 
momentarily.  Next,  the  SLR  shell  center  node  displacement  rate  increases  dramatically. 
Globally,  the  entire  displacement  field  undergoes  a  quick  reversal  defined  as  the  smqiping 
phenomenon  (2:631).  The  Donnell  shell  does  not  sn^  and  is  discussed  later  on. 
Measurement  of  Avv/Ar  fi'om  the  graph  indicates  that  the  velocity  of  the  center  node 
reaches  speeds  in  excess  of  100  miles  per  hour  when  the  shell  snqis.  Immediately  after 
the  shell  reaches  its  maximum  displacement,  the  center  node  enters  into  a  state  where  it 
vibrates  with  an  amphtude  and  period  that  decays  over  a  short  period  and  will  be  referred 
to  as  a  transient  phase  in  this  paper.  After  the  transient  phase,  the  center  node  oscillates 
about  an  equilibrium  point  which  is  near  the  static  equilibrium  load.  Kounadis  noted 
similar  displacement  characteristics  with  his  two-degree  of  freedom  pendulum  model  (8). 

Work  has  been  done  to  define  the  point  of  dynamic  instability  of  a  system.  Kounadis 
analyzed  one  and  two-degree-of-freedom  systems,  with  material  nonlinearities  (instead  of 
geometric  nonlinearities),  subjected  to  step  loading  and  defined  dynamic  buckling  as  "that 
state  for  which  an  escaped  motion  is  initiated"  (8).  The  escaped  motion  was  associated 
with  an  angular  degree  of  freedom  with  subsequent  unbounded  behavior  over  time.  The 
sudden  rise  in  velocity  as  indicated  in  Figure  4.1  is  analogous  to  the  escaped  motion 
observed  by  Kounadis  with  the  exception  that  unbounded  motion  is  not  observed  nor 
expected.  For  this  higher  degree-of-freedom  system,  Chien  and  Palazotto  (2)  defined 
the  inflection  point  on  the  displacement  curve  as  the  point  where  dynamic  instability 
initiates  for  arches  and  cylindrical  shells.  A  general  discussion  of  this  dynamic  instability 
phenomenon  from  energy  considerations  may  be  found  there.  In  short,  the  snipping 
phenomenon  is  a  sudden  transference  of  stored  potential  energy  to  kinetic  energy  at  a 
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certain  threshold  level.  This  effect  may  be  seen  at  the  early  part  of  the  displacement 
history  in  Figure  4.1.  Upon  initial  t^>piication  of  the  load,  there  is  an  increase  in  velocity. 
Soon  after,  the  average  velocity  decreases  considerably  while  tl^  peak  load  level  is 
obtained;  potential  energy  is  stored.  Continued  load  q)plication  eventually  releases  the 
stored  energy. 

It  is  interesting  to  note  that  the  8  in.  radius  shell  does  not  initiate  snapping  until 
0.0025s  after  the  maximum  load  has  been  reached  as  indicated  by  the  pairs  of  vertical 
lines  in  Figure  4.1.  This  delay  is  attributed  to  a  period  of  dynamic  instability  that 
dominates  the  behavior  even  though  the  static  coll^rse  load  level,  which  is  sufficient  for 
dynamic  collapse,  has  been  reached.  Chien  and  Palazotto  (2)  showed  that  the  onset  of 
dynamic  instability  can  be  determined  by  comparing  the  behavior  of  a  shell  with  a  load 
that  does  not  initiate  collapse  with  the  behavior  of  the  same  shell  with  a  load  that  just 
initiates  collapse.  The  onset  of  dynamic  instability  occurs  at  the  point  where  the  two 
displacement  curves  diverge.  For  the  8  in.  radius  shells,  this  divergence  takes  place  at 
approximately  t=  0.003  which  is  coincident  with  the  time  when  the  maximum  load  has 
been  achieved.  The  pre-collapsed  behavior  of  the  SLR  shell  is  similar  to  the  dotted  curve 
in  Figure  4.1.  In  general,  the  Doimell  curves  are  slightly  lower  than  the  SLR  curves  in 
the  pre-collapsed  phases).  Comparatively  speaking,  this  period  of  dynamic  instability  for 
the  8  in.  radius  shell  was  more  prolonged  than  the  shallower  arches  and  shells  studied  by 
Chien  and  Palazotto.  Therefore,  there  appears  to  be  a  relation  between  the  depth-to-chord 
{hid)  ratio,  namely  that  dynamic  instability  is  expected  to  be  more  pronounced  for  shells 
with  higher  hic  ratios.  However,  the  mechanism  for  the  delay  caused  by  the  instability 
is  not  fully  understood.  Chien  and  Palazotto  (2)  have  suggested  that  the  dynamic 
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instability  is  associated  with  a  dynamic  coupling  between  the  primary  state  fiinctitm  and 
a  secondary  state  displacement  state  known  as  parametric  resonance.  Simitses  (22)  noted 
the  potential  for  parametric  resonance  to  occur  in  a  suddenly  loaded  structure.  That  is, 
besides  the  motion  produced  by  the  forcing  function,  there  may  be  other  significant 
mo  ins  not  necessarily  in  the  direction  as  the  forcing  function.  Evidence  of  this 
phenomenon  is  presented  below. 

Figure  4.2  shows  a  sequence  of  meshes  that  depict  the  composite  shell  at  various 
times  in  the  loading  history.  It  should  be  also  noted  that  the  meshes  represent  the  datum 
surface  of  the  composite.  The  meshes  displayed  generally  correspond  to  significant  points 
on  the  center  node  displacement  curve,  i.e.,  the  inflection  point,  local  maxima  and  minima 
during  the  transient  phase,  and  the  later  stages.  The  figure  shows  that  an  initial  dimpling 
takes  place  followed  by  the  formation  of  a  deep  crease  along  the  longitudinal  direction. 
Very  soon  after,  the  shell  reaches  a  maximum  displaced  state  in  a  very  short  time,  i.e., 
the  snapping  phenomenon. 

During  the  transient  phase  (tqrproximately  t=0.0075  to  0.014  seconds).  Figure  4.2  at 
t  =  O.OOSSs  shows  that  the  free  edges  of  the  shell  undergo  secondary  sniping  motions. 
When  the  shell  bottoms  out,  the  nodal  velocities  are  minimal  and  potential  energy  is 
stored.  Since  the  center  node  is  the  only  node  that  has  direct  loading,  all  other  nodes  are 
free  to  give  up  stored  elastic  energy,  hence  the  secondary  sniping.  After  the  transients 
have  stopped,  the  displacement  field  stays  reversed  with  no  further  large  displacements 
taking  place.  Vibratory  motion  in  the  w  direction,  after  the  transient  phase  has  ceased, 
is  apparent  from  Figure  4.1,  but  there  are  also  in-plane  parametric  vibrations  as  weU.  The 
distorted  mesh  elements  at  the  various  times  in  Figure  4.2  is  visual  evidence  of  this  idea. 
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Meshes  rqpresenting  a  composite  cylindrical  shell  subjected  to  a  transverse  stq)  load  applied  to  the  center. 
SLR  fmmulation. 


To  further  clarify  this  point.  Figure  4.3  shows  a  1”  x  1"  square  element  at  various  times 
in  its  loading  histmy.  In  order  to  illustrate  the  changing  shape  of  the  element,  the  lower 
left  corners  of  the  elements  are  fixed  to  the  same  location,  i.e.,  translation  of  the  defmined 
element  with  respect  to  the  undeformed  element  is  removed.  Figure  4.3a  corresponds  to 
times  generally  during  the  transient  phase  whereas  Figure  4.3b  corresponds  to  even 
intervals  of  time  after  the  transient  phase.  Comparisons  between  the  deformed  and 
undeformed  shell  elements  give  a  qualitative  impression  of  the  strain  and  stress  fields  as 
they  change  with  time.  In  general,  one  would  surmise  that  after  collapse  the  elements 
closer  to  the  hinged  boundaries  are  primarily  in  a  state  of  tension  while  the  elements 
towards  the  middle  are  cycling  between  tensile  and  compressive  states.  A  cursory 
examination  of  the  shell  elements  also  indicates  that  distortions  are  more  evident  in  the 
V  in-plane  directions  (circumferential)  rather  than  the  u  directions.  A  check  of  the  output 
data  files  for  all  of  the  shells  studied  indicates  that  u  displacements  are  indeed  minimal 
and  may  be  due  to  rotational  effects  only.  Since  numerical  models  are  ideal  from  a 
material  and  load  applications  standpoint,  excitation  of  the  shell  in  the  longitudinal 
direction  does  not  occur  for  the  composite  layup  orientations  considered. 

The  lower  curve  in  Figure  4.1  is  a  result  from  the  Donnell  model.  The  curve 
indicates  that  the  Donnell  shell  does  not  collapse  despite  the  fact  that  the  loading  is  the 
same  as  the  SLR  shell.  Recall  that  the  collapse  load  level  was  calculated  with  a  static 
version  of  the  SLR  formulation.  This  result  indicates  that  the  Donnell  shell  has  stiffer 
effective  flexural  characteristics  than  the  higher-order  formulation.  It  is  believed  that  this 
higher  stiffness  is  a  by-product  of  the  absence  of  higher-order  terms  in  the  formulation. 
At  first  glance,  the  Donnell  shell  goes  primarily  into  a  quasiperiodic  oscillation  mode; 
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Radius  =  8  in.,  SLR  formulation.  Dotted  lines  correspond  to  the  underformed  state. 


however  this  is  n(^  the  case  as  will  be  demcmstrated  in  Section  VI.  A  display  of  the 
defcHined,  uncollapsed  mesh  is  not  provided,  but  is  similar  to  the  precollapsed  SLR  mesh 
in  Figure  4.2  at  /  »  0.0045s. 

42  Twelve  in.  Radius  Cylindrical  Shell 

Figures  4.4  and  4.S  show  the  center  node  di^lacement  charactmistics  of  a  12  in. 
radius  cylindrical  shell  using  botft  formulations.  The  curve  characteristics  of  Figures  4.4 
and  4.5  are  similar  to  the  SLR  curve  in  Figure  4.1  in  that  there  is  an  instability  inflection 
point,  a  sudden  increase  in  the  velocity  of  the  displacement  field,  a  decaying  transient 
period,  and  an  oscillation  about  an  equilibrium  level.  Unlike  the  8  in.  radius  shell,  both 
12  in.  radius  shells  collq>se  under  the  same  loading  history.  It  spears  that  either  1)  the 
increase  in  the  cylindrical  radius  has  the  effect  of  reducing  the  effective  flexural  stiffness 
(see  Appendix  A)  of  the  Donnell  shell  to  a  level  close  to  that  of  the  SLR  shell  (recall  that 
previous  results  indicate  that  the  Donnell  shell  is  stiffer  than  the  SLR  shell)  or  2)  the 
applied  load  is  well  beyond  the  levels  required  to  cause  coll^se.  In  either  case,  there  is 
a  very  small  delay  between  the  occurrence  of  the  maximum  load  and  the  onset  of 
sni4}ping.  This  result  is  further,  evidence  that  the  dynamic  instability  is  related  to  the  /i/c 
ratio  which  in  turn,  is  related  to  the  stiffness  prq)erties  of  the  shell. 

Figures  4.6  and  4.7  show  a  sequence  of  meshes  depicting  the  12  in.  radius  shell  as 
it  deforms.  The  similarities  between  the  mode  shapes  of  the  two  sets  of  figures  are  also 
{parent  and  indicate  the  good  agreement  between  the  SLR  and  Donnell  models  for  this 
shell  geometry.  The  primary  difference  between  Figures  4.6  and  4.7  is  in  the  transient 
phase  following  the  shell  collapse.  The  transient  phase  of  the  SLR  shell  has  a  decaying 
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Figure  4.4.  Center  node  displacement  of  a  composite  cylindrical  shell  subjected  to  a 
transverse  step  load  applied  to  the  center.  Radius  =  12  in.,  SLR  fcnmulation. 
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Figure  4JS.  Center  node  displacement  of  a  composite  cylindrical  shell  subjected  to  a 
transverse  step  load  applied  to  the  center.  Radius  =  12  in.,  Donnell  fmmulation. 
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re  4.6b.  Meshes  representing  a  composite  cylindrical  shell  subjected  to  a  transvose  stq>  load  applied  to  the  center, 
us  B  12  in.,  SLR  fonnulation. 
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ius  s  12  in.,  Donnell  formulation. 


amplitude  and  an  increasing  fiequency  at  the  center  node.  On  the  other  hand,  the  center 
node  of  the  Etonnell  shell  experiences  a  significant  upward  motion  at  the  end  of  the 
transient  phase.  Further,  the  entire  displacement  field  (except  at  the  boundaries) 
experiences  a  sudden  upward  movement;  this  phenomenon  is  illustrated  in  Figure  4.7  at 
t  =  0.0085s.  This  effect  is  not  apparent  in  the  SLR  models.  In  general,  the  potential  and 
kinetic  energy  distribution  of  the  Donnell  model  distributes  throughout  the  shell 
differently  than  the  SLR  model  because  energy  that  would  normally  distribute  through  all 
of  the  higher-order  nonlinear  terms  (numerically  through  solution  of  the  matrix  equations) 
is  not  possible  in  the  Donnell  model.  This  effect  shows  up  gr^hically  in  figures  for  the 
deeper  shells  in  the  following  section. 

4.3  Six  in.  Radius  Cylindrical  Shell 

For  a  deeper  6  in.  radius  cylindrical  shell,  both  model  formulations  are  observed  to 
initiate  collapse  imder  the  same  load  history  as  indicated  in  Figure  4.8.  Although  the 
effective  flexural  stifihess  of  the  Donnell  shell  is  higher  than  the  SLR  shell  (as  it  was  in 
the  8  in.  radius  shell),  the  applied  loading  is  sufficient  to  result  in  collapse  for  both  shells. 
However,  it  is  likely  that  there  is  a  small  range  of  load  levels,  for  a  given  loading  rate, 
where  only  the  SLR  shell,  and  not  the  Donnell  shell,  would  collapse.  Figure  4.8  shows 
that  the  Donnell  curve  also  experiences  a  significant  upward  displacement  near  the  end 
of  the  transient  period  which  is  a  behavior  observed  in  tlte  12  in.  radius  Donnell  shell. 
Moreover,  Figure  4.8  at  time  1=0.01 15s  is  also  similar  to  Figure  4.4  at  /=0.(X)85s  in  that 
the  displacement  field  snaps  a  small  distance  upwards.  It  is  also  apparent  fi'om  Figure 
4.8  that  the  differences  between  the  Donnell  and  SLR  formulations  are  more  conspicuous. 
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Figure  4.8.  Center  node  displacements  of  composite  cylindrical  shells  subjected  to  a 
transverse  step  load  applied  to  the  centm*.  Radius  =  6  in. 
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Prior  to  collapse,  both  formulations  predict  similar  behavior.  Also,  the  reduction  of  the 
center  nocte  displacement  rate  in  the  early  part  of  the  curve  is  more  conspicuous  when 
compared  to  shallower  shells.  Again,  collapse  occurs  a  significant  time  aftm*  the 
maximum  load  level  is  achieved  as  indicated  in  Figure  4.8.  As  Chien  and  Palazotto 
indicate  (2),  the  point  of  divergence  between  curves,  corresponding  to  a  shell  that 
eventually  collapses  and  a  shell  that  does  not  collapse,  indicates  the  initiation  of  dynamic 
instability.  A  curve  corresponding  to  a  sl^ll  that  does  not  coUqrse  was  generated  but  is 
not  shown.  As  in  the  8  in.  radius  shell,  the  dynamic  instability  initiates  at  the  same  time 
the  maximum  load  level  is  reached.  The  vertical  lir^  in  Figure  4.16  indicates  the  time 
when  the  maximum  load  level  is  reached  (the  initiation  of  dynamic  instability  possibly 
due  to  parametric  vibration). 

After  colli^se,  the  Donnell  shell  displaces  a  further  distance.  Also,  the  Donnell  shell 
has  more  oscillations  during  the  transient  phase.  Figures  4.9  and  4.10  show  the  Donnell 
and  SLR  displacement  fields  at  various  times  for  the  6  in.  radius  shell.  For  the  stages  up 
to  and  at  collapse,  the  shell  shares  are  similar.  During  the  transient  phases,  the  behaviors 
are  very  different.  In  particular,  the  Donnell  formulation  exhibits  some  peculiar  behavior. 
For  example.  Figure  4.10  shows  that  some  meshes  have  elements  along  the  boundaries 
that  have  stretched  to  several  times  its  original  size,  while  the  elements  along  the 
centerline  have  compressed  to  very  small  sizes.  In  particular,  the  Dormell  mesh  at  r  = 
0.(X)95s  takes  on  an  odd  t^rpearance.  This  behavior  shows,  with  more  clarity  (from  a 
graphical  viewpoint),  that  the  distribution  of  energy  in  the  Donnell  model  is  significantly 
differently  than  the  SLR  model.  The  Donnell  solution  has  the  elements  along  the 
boundary  in  a  higher  state  of  tension  (indicated  by  stretched  elements)  and  the  elements 
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ire  4.9b.  Meshes  representing  a  composite  cylindrical  shell  subjected  to  a  transverse  st^  load  applied  to  the  center, 
ius  s  6  in.,  SLR  formulation. 
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along  the  center  line  in  an  higher  state  of  compression  (indicated  by  compressed 
elements),  when  compared  to  the  SLR  model  during  the  transient  phase.  The  elements 
in  the  SLR  model  tend  to  be  more  evenly  proportioned  relative  to  the  original  undeftnmed 
element  shapes.  After  the  transient  period,  both  models  show  that  the  elements  are 
distributed  more  evenly  in  size. 

4.4  Four  in.  Radius  Cylindrical  Shell 

Figure  4.11  shows  the  center  node  displacement  curves  of  4  in.  radius  cylindrical 
shells  using  the  SLR  and  Donnell  formulations.  The  disagreement  between  both  models 
is  even  more  ^parent  after  collapse,  although  the  curve  characteristics  as  discussed  in  the 
preceding  paragrt^hs  are  the  same.  The  inflection  point,  which  is  an  indication  of  the 
onset  of  instability,  is  not  as  distinguishable  as  it  is  in  the  shallower  shells;  this  may  be 
due  to  the  maximum  load  level  occurring  in  a  more  gradual  fashion  as  indicated  in  the 
loading  curve.  The  SLR  shell  sn^s  sooner  than  the  Donnell  shell  since  the  higher 
number  of  rotational  terms  gives  the  SLR  shell  more  flexibility.  Figures  4.12  and  4.13 
display  the  deformed  geometries  of  the  Donnell  and  SLR  shells.  Despite  the  fact  that  the 
models  assume  perfectly  elastic  behavior  with  no  material  failure,  the  odd  post-coll^sed 
meshes  suggest  that  the  SLR  and  Donnell  formulations  may  only  be  applicable  for  the 
analysis  of  deep  cylindrical  shells  in  the  pre-collapsed  stages.  The  post-collapsed 
geometries  simply  do  not  appear  to  be  consistent  with  physical  expectations.  A 
comparison  between  deformed  shell  elements  with  the  undeformed  elements  suggests  that 
inordinate  levels  of  strain  are  taking  place  in  both  models,  especially  near  the  hinged 
boundaries.  In  addition,  ^parent  graphical  aberrations  in  the  meshes,  such  as  the  ones 
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Figure  4. 11.  Center  node  displacements  of  composite  cylindrical  shells  subjected  to 
a  transverse  step  load  applied  to  the  center.  Radius  =  4  in. 
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Figure  4.12.  Meshes  representing  a  composite  cylindrical  shell  subjected  to  a  transverse  step  load  applied  to  the  center. 
Radius  =  4  in.,  SLR  formulation. 
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Radius  =  4  in.,  Donnell  formulation. 


displayed  at  /  =  0.0092s  in  Figure  4.12  and  three  of  the  meshes  in  Figure  4.13  are  a  result 
of  portions  of  the  laminate  passing  through  its  symmetric  half.  Further,  it  {q>pears  that 
there  are  severe  sl(^  discontinuities  between  many  of  the  elements'.  Still,  compatibility 
and  potential  energy  constraints  are  satistied  even  though  the  shiq)es  do  not  correlate  well 
with  a  real  system.  This  problem  indicates  a  limitation  of  the  Donnell  and  SLR 
formulations  at  these  levels  of  rotations.  It  is  possible  that  refining  the  mesh  may 
improve  the  results.  The  initiation  of  somewhat  similar  trends  in  the  6  in.  radius  shell 
may  also  appear  to  question  some  of  the  peculiar  features  observed  in  that  solution. 

45  Twenty  in.  Radius  Cylindrical  Shell 

For  cylindrical  shells  shallower  than  the  12  in  radius  shell,  the  SLR  and  Donnell 
formulations  produce  very  similar  results.  Figures  4.14  and  4.15  show  that  the  center 
node  displacements  for  a  20  in.  radius  cylindrical  shell  have  very  similar  general 
characteristics.  The  differences  between  the  two  curves  would  only  be  one  of  precise 
detail.  The  displacement  curve  has  many  of  the  characteristics  noted  earlier,  i.e.,  an 
inflection  before  coll^^se,  a  snapping  action,  a  transient  period,  and  an  oscillation  near 
the  static  equilibrium  solution.  As  in  all  previous  cases,  the  initiation  of  collapse  occurs 
slightly  after  the  peak  load  is  attained,  but  this  period  of  dynamic  instability  is  very  short. 

A  comparison  of  the  deformed  meshes  between  both  models  indicates  no  significant 
differences  between  the  two  formulations.  The  coU^sed  shell  of  Figure  4.16  could 
represent  the  coll^sed  sht^  of  either  the  Donnell  or  SLR  models.  As  in  L  ^  other  cases, 

^It  should  be  noted  that  this  discontinuity  is  primarily  a  limitation  in  the  software  package 
used  to  display  the  meshes.  The  node  defining  an  element  are  connected  by  straight  liiKS  only. 
Therefore  interelement  slope  compatibility  is  not  accurately  rqtresented. 
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Figure  4.14.  Center  node  displacement  of  a  composite  cyUndrical  shell  subjected  to  a 
transverse  step  load  q)plied  to  the  center.  Radius  =  12,  SLR  formulation. 


Figure  4.15.  Center  node  displacement  of  a  composite  cyhndrical  shell  subjected  to  a 
transverse  step  load  applied  to  the  center.  Radius  =  20,  Donnell  formulation. 
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one  may  observe  slight  deformations  of  the  shell  elements  indicating  the  presence  of  in¬ 
plane  vibrations.  However,  these  changes  are  very  slight  from  a  graphical  standpoint,  and 
for  this  reason,  no  other  collapsed  meshes  are  displayed. 

4.6  Sixty  in.  Radius  CylincU-ical  Shell 

For  a  60  in.  radius  shell,  there  appears  to  be  no  discemable  deflection  point;  the  shell 
appears  to  collapse  immediately  as  indicated  in  the  center  node  displacement  history  of 
Figure  4.17.  One  significant  difference  between  this  displacement  curve  and  the  previous 
curves  is  an  absence  of  very  high  firequency  oscillations  superimposed  on  any  portion  of 
the  waveform.  Due  to  the  higher  radius  of  this  shell,  when  compared  to  the  smaller  radii 
shells,  the  behavior  is  more  oriented  to  plate  response  as  will  be  subsequently  discussed. 
This  figure  also  shows  that  the  SLR  and  Donnell  models  are  in  good  agreement  during 
the  first  one-and-a  half  cycles  but  that  the  difierent  flexural  behaviors  predicted  by  the 
two  formulations  result  in  a  slight  phase  difference  at  latter  cycles. 

4.7  Flat  Plate 

A  cylindrical  shell  with  an  arbitrarily  large  radius  of  900,000  in.  is  assumed  to  behave 
as  a  fiat  plate.  For  a  flat  plate,  there  is  no  flexural  stiffiiess  for  a  transverse  load  to 
overcome  and  hence  no  snapping.  The  stress  field  is  predominantly  in  tension  upon 
immediate  application  of  the  load.  Various  excitation  loads  were  applied  to  the  flat  plate. 
Figure  4.18  shows  loading  histories  applied  to  a  flat  plate  where  the  maximum  load  level 
is  the  same  as  the  load  applied  to  the  20  in.  radius  shell,  but  with  loading  rates  that  differ. 
The  motivation  was  to  determine  the  effect  of  loading  rate  on  plate  response.  Figures 
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0.0  seconds  collapsed 


Figure  4.16.  Meshes  depicting  a  composite  cylindrical  shell  in  an  unloaded  and  coll^sed 
state.  Radius  =  20  in. 


Figure  4.17.  Center  node  displacements  of  composite  cylindrical  shell  subjected  to  a 
transverse  point  load  applied  to  the  cento".  Radius  =  60  in. 
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4.19,  4.20,  4.21,  and  4.22  show  the  center  node  displacement  of  plates  that  result  from 
these  loads.  The  plate  subjected  to  a  higher  loading  rate  undergoes  larger  deflections  and 
contains  higher  frequency  components  wt^n  compared  to  the  plates  with  the  lower 
loading  rates  (this  is  further  shown  in  Section  6.).  However,  tl^  higher  frequency 
components  are  at  a  much  lower  frequency  wlwn  compared  to  the  cylindrical  shells  with 
radii  less  than  60  in.  (this  is  shown  in  Section  6  also).  The  plates  with  the  lower  loading 
rates  also  t^jpear  to  exhibit  hannonic  or  quasiperiodic  motion  after  the  tq^plied  load 
reaches  its  maximum  level;  however,  this  is  not  the  case  as  will  be  demonstrated  in 
Section  VI.  Figure  4.23  is  another  example  which  shows,  with  more  clarity,  the  effect 
that  initial  loading  rate  has  on  the  behavior  of  the  composite  plate.  This  figure  represents 
the  load  ^plied  to  the  60  in.  radius  cylindrical  shell  with  the  exception  that  the  various 
loading  rates  are  much  lower  than  any  of  tte  rates  used  thus  far. 

Figure  4.24  shows  the  center  node  displacement  of  a  flat  plate  subjected  to  the 
loading  history  of  the  8  in.  radius  sheU  (see  Figure  3.2).  Higher  deflections,  larger 
oscillation  amplitudes,  and  an  increase  in  high  frequency  components  are  apparent  when 
compared  to  the  previous  plate  examples.  Compared  to  the  8  in.  radius  shell,  the  plate 
does  not  exhibit  similar  high  frequency  behavior  (see  Figure  4.1).  This  effect  may  be 
due  the  absence  of  the  snapping  phenomenon  and  the  constrained  motion  of  a  vibrating 
plate.  The  snapping  phenomenon  is  associated  with  the  high  inotia  values  due  to  the 
sudden  change  in  the  displacement  field.  Moreover,  after  collapse  the  cylindrical  shell 
elements  vibrate  with  in-plane  and  out-of-plane  motions.  The  plate  elements,  on  the  other 
hand,  probably  vibrate  primarily  in  the  vertical  directions.  Therefore,  the  high  frequency 
character  of  the  collapsed  shell  center  node  may  be  due  in  part  to  the  disturbance  caused 
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Figure  4.18.  Various  loading  curves  applied  to  flat  composite  plates. 


Time,  seconds 


Figure  4.19.  Center  node  displacement  of  a  composite  plate  subjected  to  a  transvme 
step  load  applied  to  the  center.  Initial  loading  rate  of  3.16(10^)  Ibs/s 
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Figure  421.  Center  node  displacement  of  a  composite  plate  subjected  to  a  transverse 
step  load  applied  to  the  center.  Initial  loading  rate  of  1.58(10^)  Ibs/s 
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Figure  4JI2.  Center  node  displacement  of  a  composite  plate  subjected  to  a  transverse 
step  load  applied  to  the  center.  Initial  loading  rate  of  7.90(10*)  Ibs/s 


Figure  423.  Effect  of  loading  rate  on  the  center  node  displacement  of  a  flat  composite 
plate  subjected  to  a  transverse  step  load  iqjplied  to  the  center.  Max  load  =  352  lbs. 
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Figure  4.24.  Center  node  displacement  of  a  flat  composite  plate  subjected  to  a 
transverse  step  load  applied  to  the  center.  Initial  loat^g  rate  of  1.4(l(f)  Ibs/s,  max 
load  of  4232  lbs. 
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by  the  other  elements.  However,  no  data  is  available  to  suf^rt  this  su{^)osition  and 
follow-up  work  is  indicated. 

Figure  4.25  shows  exaggerated  views  of  some  of  the  displacement  fields 
corresponding  to  the  plate  with  the  heavier  loading  (a  correct  scaling  produced  meshes 
with  features  that  were  difficult  to  distinguish).  The  tigure  indicates  that  the  free  edges 
of  the  flat  plate  undergo  significant  motions  as  in  the  deep  shell  cases,  but  these  effects 
are  diminished  in  intensity.  Thus,  as  the  shell  curvature  is  reduced,  the  edge  effects 
restrains  the  free  edge  motion  in  relation  to  the  displacement  at  the  center  node. 

An  unexpected  outcome  of  the  plate  analysis  is  that  there  was  no  significant 
difference  between  the  SLR  and  von  Karman  (i.e.,  a  special  case  of  the  Donnell 
formulation  where  all  radius  values  are  arbitrarily  large)  results  for  a  flat  plate. 
Comparison  between  the  output  files  of  the  two  formulations  did  reveal  very  minute 
numerical  differences,  but  not  enough  to  afiect  the  displayed  waveforms.  Thus,  all  of  the 
previous  plate  results  are  {q)plicable  to  both  the  SLR  and  von  Karman  formulations.  This 
result  implies  that  the  higher-order  nonlinear  terms  are  effectively  removed  when  the 
value  of  Rj  is  arbitrarily  large.  A  check  of  the  strain  equations  (14:44)  indicates  that  Rj 
q)pears  in  the  denominators  of  most  (if  not  all)  of  the  higher-order  nonlinear  terms 
effectively  removing  them  from  the  equations. 
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Figure  4.25b.  Meshes  i^resenting  a  flat  composite  plate  subjected  to  a  transverse  step  load.  Vertical  displacements  scaled 
up  by  a  factor  of  4. 


V  NUMERICAL  STABIUTY 


The  time-step  selection  scheme  outhned  in  the  Section  in  results  in  numerically  stable 
output  for  all  of  the  Donnell  models  examined.  On  the  other  hand,  for  SLR  cylindrical 
shells  with  radii  less  than  8  inches,  the  recommended  post-collapse  time  step  will 
eventually  result  in  the  accumulation  of  numerical  errors.  It  is  believed  that  since  the 
SLR  formulation  contains  far  more  terms  than  the  Donnell  formulation,  the  possibility  of 
propagating  eiJ’Oi's  becomes  higher.  Figures  5.1  and  5.2  show  where  numerical 
instabilities  occur  for  a  given  time  step.  Despite  the  fact  that  an  erratic  solution  is 
predicted,  which  implies  numerical  convergence,  the  solution  is  not  assumed  to  be  correct. 
Accuracy  of  this  unstable  portion  can  be  checked  by  reducing  the  time  step  further.  As 
the  figures  indicate,  a  slightly  smaller  time  step  results  in  completely  different 
displacement  values  (upper  portions  of  the  curves)  when  compared  to  the  numerical 
instabilities.  Further  time  step  reductions  do  not  result  in  ^preciable  changes  in  the 
upper  portions  of  the  curves,  so  the  solution  is  likely  to  be  correct.  The  figures  suggest 
that  periodically  reducing  the  time  step  will  result  in  numerically  stable  ou^ut.  However, 
the  drawback  to  this  procedure  is  in  the  significant  computational  times  required.  For 
example,  a  fully  dedicated  SPARC  2  Sun  system  required  12  days  to  compute  that  portion 
of  the  SLR  curve  in  Figure  5.2  where  a  time  step  of  2.0(10'®)  seconds  was  required  for 
stable  results.  Moreover,  it  is  not  clear  when  a  change  to  the  smaller  time  step  should 
be  initiated.  An  alternative  would  be  to  require  improved  accuracy  by  lowering  the 
convergence  tolerance  requirements,  but  it  is  also  not  clear  if  this  procedure  would  reduce 
or  increase  the  computational  time.  Although  it  is  recognized  that  the  SLR  and  Donnell 
formulations  are  not  applicable  to  modelling  the  post-collq)sed  behavior  of  the  deeper 
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Figure  5.1  Numerical  instability  in  tlK  center  node  displacement  of  a  composite 
cylindrical  shell  subjected  to  a  transverse  step  load  applied  to  the  center.  Radius  =  6  in., 
SLR  formulation. 


Figure  5.2  Numerical  instability  in  the  center  node  displacement  of  a  composite 
cylindrical  shell  subjected  to  a  transverse  step  load  appUed  to  the  center.  Radius  =  4  in., 
SLR  formulation. 
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shell  geometries  for  the  meshes  selected  (it  is  not  clear  if  changing  the  mesh  would 
produce  better  results),  the  findings  do  indicate  that  the  Newton-Raphson  and  Beta-m 
methods  will  produce  numerically  stable  results  for  deep  shells  as  long  as  the  time  step 
is  sufficiently  small.  Therefore,  these  iterative  methods  should  be  useful  in  producing 
stable  results  for  formulations  incorporating  large  deflection  and  large  rotation  features. 

Numerical  instability  is  also  evident  from  an  examination  of  the  corresponding  mesh. 
Figure  5.3  shows  an  example  of  a  numerically  unstable  mesh  solution.  Compatibility  and 
potential  energy  constraints  are  still  satisfied,  but  the  accumulation  of  numerical  errors 
results  in  a  totally  unrealistic  displacement  field. 


Figure  5.3  Numerically  unstable  solution  of  a  deep  composite  cylindrical  shell  subjected 
to  a  transverse  point  load  applied  to  the  center. 
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VI  CHAOTIC  BEHAVIOR 


There  is  ample  evidence  to  show  that  undamped  cylii^cal  shells  exhibit  chaotic 
behavior  after  collapse.  Figure  6.1  shows  the  phase-plane  trajectory  of  the  collapsed 
response  of  the  20  in.  radius  cylindrical  shell.  This  graph  indicates  that  the  orbits  never 
close  or  repeat  and  that  the  trajectory  fills  up  a  portion  of  the  phase  space  which  are 
characteiisitic  of  chaotic  motion.  Figure  6.2  shows  a  Poincar6  map  of  the  same  shell 
since  the  applied  force  was  not  of  the  periodic  type,  the  sampling  period  T  was  chosen 
to  be  an  average  of  the  time  between  dominant  amplitude  peaks.  Figure  6.3  shows  an 
expanded  portion  of  the  SLR  displacement  curve  (and  the  data  point  resolution  used  in 
plotting  the  curve).  The  highly  dynamic  changes  of  the  curve  illustrates  the  dilemma  in 
choosing  an  appropriate  sampling  time;  if  a  periodic  forcing  function  was  ^phed,  the 
sampling  time  could  have  been  chosen  by  the  criteria  discussed  in  Section  2.3.  Figure 
6.2  indicates  that  the  a  cloud  of  unorganized  points  inhabits  a  portion  of  the  phase  space 
plane;  this  is  indicative  of  the  existence  of  a  strange  attractor  (as  defined  in  Section  2.3) 
which  draws  the  trajectories  to  a  region  in  the  phase  space,  but  due  to  the  lack  of 
dissipation,  fractal  structures  (as  defined  in  Section  2.3)  are  not  expected. 

Figure  6.4  shows  the  phase-plane  trajectory  of  a  flat  plate  subjected  to  the  same 
loading  history  applied  to  the  20  in.  radius  cylindrical  shell.  The  orbits  have  more 
ordered  trajectories  when  compared  to  Figure  6.1.  Since  the  curves  are  not  closed  and 
repetitive,  the  motion  has  the  potential  to  be  classified  as  chaotic.  Moreover,  a  Poincard 
map  may  indicate  the  presence  of  fi'actal-like  structures  since  the  trajectories  tend  to 
congregate  about  similar  orbital  paths. 

Unfortunately,  the  present  computational  resources  limits  the  amount  of  available  data 
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Figure  6.1.  Phase-plane  trajectory  for  chaotic  motion  of  a  coU^sed  20  in.  radius 
composite  cylindrical  shell  subjected  to  a  transverse  point  load  applied  to  the  center. 
Shell  is  undamped. 


Figure  6.2.  Poincard  map  of  chaotic  motion  of  a  collapsed  20  in.  radius  composite 
cylindrical  shell  with  no  damping. 
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Figure  6.3.  Portion  of  the  center  node  displacement  history  of  a  collapsed  composite 
cylindrical  shell  subjected  to  a  transverse  point  load.  Radius  =  20  in.,  SLR  formulation. 


Figure  6.4.  Phase-plane  trajectory  for  chaotic  motion  of  a  flat  composite  plate  subjected 
to  a  transverse  point  load  tqjphed  to  the  center. 


that  may  be  used  to  further  clarify  the  existence  of  chaotic  phenomenon.  For  example, 
since  one  point  on  a  Poincar6  map  represents  a  minimum  of  one  cycle,  hundreds  or 
thousands  of  cycles  are  needed  to  adequately  capture  the  character  of  a  strange  attractor. 
As  Figure  4.20  indicates,  there  are  at  most  1 1  cycles.  Nearly  3  days  of  computational 
time  were  required  to  generate  the  data  for  that  figure.  Further,  it  would  be  desirable  to 
quantify  existence  of  chaotic  behavior  with  the  calculation  of  Lyapunov  exponents,  but 
the  time  required  to  generate  the  necessary  amount  of  data  was  not  available  towards  the 
conclusion  of  this  work.  Zeng,  Eykholt,  and  Pielke  (29)  have  demonstrated  that  a 
minimum  of  about  5(XX)  data  points  in  three-  or  four-dimensional  phase  space,  without 
transients,  is  needed  to  obtain  a  reasonable  Lyapunov  value.  The  largest  number  of  data 
points  generated  for  one  particular  shell  geometry  in  this  work  numbered  no  more  than 
2500  data  sets  (including  transients)  and  required  in  excess  of  2  weeks.  Attempts  were 
made  to  utilize  Wolfs  algorithm  (27)  to  determine  the  largest  Lyq)unov  exponent  for  the 
collapsed  12  and  20  in  radius  shells  since  there  appeared  to  be  several  hundred  cycles 
available.  However,  Wolfs  algorithm  is  generally  oriented  to  systems  with  dissipation. 
The  highly  dynamic  behavior  of  the  data  for  the  two  shell  geometries  did  not  lead  to  a 
convergence  to  a  Ly{q)unov  exponent  The  calculation  of  Lyapunov  exponents  is  based 
on  exponential  growth  between  a  reference  trajectory  (called  a  fiduciary  by  Wolf  et  al. 
(27)  aiKl  a  point  on  a  nearby  trajectory.  When  the  growth  between  the  two  trajectories 
grows  at  a  rate  which  is  larger  than  an  exponential  rate,  the  algorithm  seeks  out  a  new 
nearby  trajectory.  As  Figure  6.3  might  suggest,  the  erratic  behavior  of  the  signal  may  not 
be  suitable  for  the  algorithm  since  large  changes  take  place  ov^  a  very  short  period. 

A  cursory  examination  of  the  Donnell  curve  in  Figure  4. 1  would  suggest  that  the  pre- 
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coll{^)sed  shell  center-node  displacement  histcny  is  a  result  of  the  supeiposition  of 
harmonic  motions.  However,  a  truly  quasipehodic  signal  would  have  closed  and  exactly 
repeating  orbits.  Figure  6.5  is  a  phase-plane  trajectory  of  the  pre-colliq>sed  8  in.  radius 
shell  and  shows  harmonic  or  quasipehodic  motion  is  not  present.  Chien  and  Palazotto  (2) 
noted  harmonic  motion  of  arches  subjected  to  step  loads.  However,  tte  applied  loads 
were  significantly  less  than  the  collapse  load.  Greer  (6)  also  noted  harmonic  motion  in 
composite  cylindrical  shells  when  the  applied  load  was  far  less  than  the  collapse  load. 
He  concluded  that  the  nonlinear  portions  of  the  SLR  formulation  played  a  minimal  part 
in  the  displacement  solution  for  that  case.  Yet  at  the  same  time,  he  noted  that  the  motion 
of  a  free  edge  node  exhibited  chaotic  tendencies.  Figure  6.5  suggests  that  chaotic 
tendencies  are  present  in  the  8  in.  radius  uncoll^sed  shell.  A  possible  explanation  is 
that  since  the  applied  load  is  close  to  the  collapse  load,  nonlinear  terms  are  significant 
enough  to  influence  the  solutions  and  hence  throw  the  trajectory  off  a  harmonic  path. 

Figure  6.6  is  a  pseudo-phase-plane  trajectory  of  the  pre-collapsed  behavior  of  an  8 
in.  radius  Donnell  shell.  The  similarities  between  this  figure  and  the  previous  figure 
further  supports  the  possibility  of  chaotic  behavior  since  the  orbits  are  nonrepetitive.  The 
similarities  between  the  two  figures  also  validates  that  a  single  time  dq)endent  variable 
may  provide  enough  information  to  classify  the  existence  of  chaotic  or  non-chaotic 
behavior,  especially  when  only  one  vaiiable  is  measured  (numerically  or  experimentally). 
Moon  (10:43)  suggests  that  one  may  numerically  integrate  (if  velocity,  for  example,  is  the 
only  measured  parameter)  or  differentiate  (if  displacement,  for  example,  is  the  only 
measured  parameter)  in  order  to  obtain  the  necessary  information  to  plot  a  phase-plane 
history  or  a  Poincar6  mq).  Again,  due  to  the  tendency  of  the  trajectories  to  ccmgregate 
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Figure  6JS.  Hiase-plane  trajectory  of  an  uiKoUapsed  composite  cylindrical  shell  subjected 
to  a  transverse  step  load  applied  to  the  center.  Radius  =  8  in.,  Donnell  formulation. 


Figure  6.6.  Pseudo-phase-plane  trajectory  for  an  uncollapsed  composite  cylindrical  shell 
subjected  to  a  transverse  step  load  applied  to  the  center.  Ra^us  =  8  in.,  Donnell 
formulation. 
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about  similar  orbital  paths,  a  PoincarS  map  might  reveal  the  presence  of  fractal  structures 
if  enough  data  points  were  available. 

A  spectrum  analysis  of  the  displacement  history  suggests  the  presence  of  chaotic 
behavior.  Figures  6.7  and  6.8  show  the  frequency  spectrums  of  the  12  in.  radius  shells 
after  colli4>se.  The  response  over  a  broad  range  of  frequencies  implies  chaotic  behavior. 
The  results  also  show  that  the  SLR  and  Donnell  formulations  yields  similar  coUapsed- 
shell  frequency  responses.  In  both  cases,  tl^re  is  a  gradual  attenuation  of  signal  strength 
as  frequency  increases.  Figure  6.9  shows  the  frequency  spectrum  of  the  20  in.  radius 
shell  using  the  SLR  formulation.  The  Donnell  model  yielded  similar  results  and  is  not 
presented.  Figure  6.10  is  the  frequency  spectrum  of  the  uncoU^sed  8  in.  radius 
cylindrical  shell.  In  a  periodic  or  quasiperiodic  signal,  there  would  be  discrete  peaks 
corresponding  to  the  signal  frequencies  separated  by  flat  responses  in  between.  The 
figure  shows  that  the  predominant  signals  form  at  the  lower  end  of  Uk  scale  as  expected 
since  the  displacement  waveform  is  composed  of  lower  frequency  signals.  Instead  of  a 
flat  response  after  the  low  frequencies,  the  signal  strength  gradually  tapers  off,  i.e.,  a 
weak,  but  broad  spectrum  of  frequencies  nonetheless.  This  characteristic  may  be  further 
evidence  that  a  shell  may  display  chaotic  behavior  without  collapsing.  A  frequency 
analysis  of  the  8  in.  radius  collapsed  shell  using  the  SLR  formulation  was  not  carried  out 
due  to  the  limited  amount  of  data  available  corresponding  to  the  collapsed  state. 

Figures  6.11,  6.12,  6.13,  and  6.14  show  the  phase-plane  plots  of  three  of  the 
composite  plates  subjected  to  step  loads  that  have  different  inial  loading  rates.  Figure 
6.14  is  the  same  as  Figure  6.13  with  the  exception  that  the  scale  is  expanded  to  reveal 
more  details  of  the  phase-plane  trajectories.  The  figures  imply  that  chaotic  motion  is 
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Figure  6.7.  Frequency  spectnun  of  the  center  node  of  a  coUapi^ed  composite  cylindrical 
shell  subjected  to  a  transverse  point  load  apphed  to  the  center  node.  Radius  =  12  in.,  SLR 
formulation. 


Figure  6.8.  Frequency  spectrum  of  the  center  node  of  a  collapsed  composite  cylindrical 
shell  subjected  to  a  transverse  step  load  applied  to  the  center.  Radius=12  in.,  Donnell 
formulation. 
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Figure  6.9  Frequency  spectrum  of  the  center  node  of  a  coU^sed  composite  cylindrical 
shell  subjected  to  a  transverse  step  load  applied  to  the  center.  Radius  =  20  in.,  SLR 
formulation. 


Figure  6.10  Frequency  spectrum  of  the  center  node  of  an  uncollq^sed  composite 
cylindrical  shell  subjected  to  a  transverse  point  load  fq)plied  to  the  center.  Radius  =  8  in., 
Donnell  formulation. 
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Figure  6.11.  Phase-plane  trajectory  for  chaotic  motion  of  a  flat  composite  plate 
subjected  to  a  transverse  step  load  implied  to  die  center.  Initial  loading  rate  of  3.16(10^) 
Ibs/s,  max  load  of  1,592  lbs. 


Figure  6.12  Phase-plane  trajectory  for  chaotic  motion  of  a  flat  composite  plate  subjected 
to  a  transverse  step  load  applied  to  the  center.  Initial  loading  rate  of  7.9(10*)  Ibs/s,  max 
load  of  1,592  lbs.. 
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Figure  6.13.  Phase-plane  trajectory  for  chaotic  motion  of  a  flat  composite  plate 
subjected  to  a  transverse  step  load  ^phed  to  the  center.  Initial  loading  rate  of 
7.9(10^)  Ibs/s,  max  load  of  1.592  lbs. 


Figure  6.14.  Scaled  up  version  of  Figure  6.13  revealing  similar  orbital  trajectories. 
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present.  They  also  demonstrate  that  diere  is  a  sensitive  dq^endence  on  initial  conditions, 
i.e.,  8  definition  of  chaos.  Moreover,  as  the  loading  rate  decreases,  the  general  path  of 
all  the  trajectories  becomes  closer.  This  phenomenon  suggests  that  if  the  load  is  i4)plied 
at  a  much  slower  rate,  the  phase-plane  orbit  would  tend  to  a  point  corresponding  to  the 
static  equilibrium  state,  i.e.,  a  quasistatic  application  of  the  load.  Figures  6.15,  6.16,  and 
6.17  show  the  ftequency  responses  of  the  flat  plate  center  nodes  and  suggest  chaotic 
behavior  is  present.  As  alluded  to  earlier,  the  flat  plate  center  node  oscillates  at  generally 
lower  frequencies  when  compared  against  the  cylindrical  shell  behaviors. 

Since  the  phase-plane  orbits  of  the  plates  appeared  to  have  more  ordered  trajectories, 
it  may  be  possible  that  a  strange  attractor  is  at  woiic.  Figure  6.18  is  a  Poincar6  map  of 
the  plate  loaded  with  the  lowest  loading  rate  (see  Figures  4.22  and  6.14  ).  A  strange 
attractor  appears  to  be  influencing  the  orbits  since  the  Poincar6  points  are  congregated  in 
a  region  in  the  phase-plane,  but  the  limited  number  of  data  points  does  not  allow  the 
identifrcation  of  fractai-like  structures.  If  fractal  patterns  were  revealed,  the  predominant 
state  of  plate  tension,  instead  of  damping,  might  be  the  mechanism  by  which  this  occurs. 
Figure  6.19  is  the  displacement  curve  used  to  plot  the  Poincar6  map.  This  data  was  also 
used  in  an  attempt  to  generate  Lyapunov  exponents.  Wolfs  (27)  algorithm  was  used. 
Input  parameters  to  the  algorithm  were  based  on  some  of  his  suggested  guidelines.  For 
example,  he  suggested  that  an  embedding  dimension  m  (the  dimension  where  the  attractor 
exists)  of  3  is  often  sufficient  and  that  the  product  rm  be  approximately  equal  to  the 
mean  orbital  period.  Figure  6.20  shows  the  results.  The  two  curves  correspond  to  rm 
values  that  are  slightly  more  and  slightly  less  than  the  mean  orbital  period;  a  curve 
corresponding  to  a  value  of  rm  equal  to  the  mean  orbital  period  was  not  included  due  to 
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Figure  6.15.  Frequency  spectrum  of  the  center  node  of  a  flat  composite  plate  subjected 
to  a  transvm^  step  load  applied  to  the  center.  Initial  loading  rate  of  3.16(10^),  max.  load 
of  1.592  lbs. 


Figure  6.16.  Frequency  spectrum  of  the  center  node  of  a  flat  composite  plate  subjected 
to  a  transverse  step  load  applied  to  the  center.  Initial  loading  rate  of  7.9(10^)  Ibs/s  max. 
load  of  1,592  lbs. 
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Figure  6.17.  Frequency  spectrum  of  tlK  center  node  of  a  flat  composite  plate  subjected 
to  a  transverse  step  load  applied  to  the  center.  Initial  loading  rate  of  7.9(10*)  Ibs/s,  max. 
load  of  U92  lbs. 


Figure  6.18.  Poincare  mt^  of  a  flat  composite  plate  subjected  to  a  transverse  point  load 
applied  to  the  center.  Initial  loading  rate  of  7.9(10*)  Ibs/s,  max.  load  of  1.592  lbs. 


Figure  6.19.  Center  node  deflection  of  a  flat  composite  plate  subjected  to  a  transverse 
step  load  applied  to  the  center.  Initial  loading  rate  of  7.90(10*)  Ibs/s,  max.  load  of  1,592 


lbs. 


Figure  6.20.  Possible  convergence  of  the  Lyapunov  exponent  of  a  transversely  loaded 
flat  composite  plate  (see  Figure  6.19). 


the  confusing  overlapping  curves  that  resulted.  The  two  curves  iQq)ear  to  oscillate  about 
two  different  positive  steady  state  values.  Ly^unov  exponents  measure  the  rate  at  which 
system  processes  create  or  destroy  information  and  may  be  expressed  as  bits  per  second 
(27).  The  value  of  bits/s  refers  to  the  predictive  accuracy  of  die  data.  For  example,  if 
the  Lyapunov  exponent  was  actually  2.(XX)  bits/s  and  if  an  initial  point  were  known  with 
an  accuracy  of  one  part  per  million  (20  bits),  future  behavior  cannot  be  predicted  after 
(20  bits)/(2,(XX)  bits/s)  =  0.01  seconds,  which  corresponds  to  approximately  S  orbits 
(qiproximately  48  data  points  equally  spaced  by  4.(X10^)s  in  this  case).  The  actual 
accuracy  of  the  data  in  parts  per  million  is  not  known  at  this  time.  Nevmiheless,  the  two 
curves  show  that  the  Ly^unov  exponents  are  expected  to  remain  positive.  If  so,  this 
quantitatively  defines  the  plate  behavior  as  chaotic.  Longer  computer  runs  would  be 
needed  to  assure  that  the  exponents  do  not  suddenly  become  negative. 
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VII  SUMMARY.  CONCLUSIONS.  AND  RECOMMENDATIONS 


7.1  Summary 

The  (xnnposite  plate  and  cylindrical  shell  models  examined  in  this  research  were 
assumed  to  remain  perfectly  elastic  widi  no  dissipati<Mi  and  no  material  failure.  S  distinct 
stages  were  noted  in  the  displacement  history  of  the  cylindrical  shells  subjected  to  a 
tramverse  st^  load  ai^lied  to  the  center.  First  there  is  the  initial  increase  in  tiie  center 
node  velocity  resulting  in  a  dimpling  of  the  shell.  Second,  there  is  a  reduction  in  the 
displacement  rate  which  may  or  may  not  be  very  distinct  dq)ending  on  the  shell  radius. 
This  stage  is  manifested  in  the  creation  of  a  trough  or  crease  along  the  shell  centerline. 
Third,  shell  collapse  occurs  which  is  characterized  by  a  sudden  velocity  increase  of  the 
entire  displacement  field  (with  the  excq>tion  of  the  hinged  boundaries).  Fourth,  there  is 
a  period  of  highly  dynamic  motion  that  d^ys  witii  time.  Finally,  the  shell  settles  into 
a  stage  where  relatively  large  motions  no  longer  occur.  The  shell  is  near  the  static 
equilibrium  solution,  but  low  amplitude,  high  frequency  vibrations  are  taking  place  both 
in-plane  and  out  of  plane.  Plates  only  exhibited  the  last  two  behaviors  after  maximum 
centm*  node  deflection  occurred. 

One  characteristic  briefly  noted  in  earlier  work  with  the  DSHELL  code  is  the  large 
secondary  motions  that  are  associated  witii  the  short  but  highly  dynamic  phase  after  shell 
collapse  in  response  to  the  transversely  applied  point  load.  This  behavior  was  primarily 
noted  in  both  SLR  and  Donnell  formulations  for  most  of  the  shells.  Shallower  shells  did 
not  exhibit  tins  feature  because  lighter  loadings  did  not  provide  enough  energy  to  the 
systm  to  evoke  tins  type  of  response.  Plates  subjected  to  higher  loadings  exhibited  tins 
response  but  the  effects  were  greatly  diminished. 
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Graphical  mesh  representations  can  reveal  impmtant  qualitative  information  about 
plates  and  shells  at  a  glance.  The  most  obvious  information  is  the  overall  geometric 
cmifigurations  that  occur  when  the  load  is  apphed.  Another  feature  is  a  general 
impression  of  the  stress/strain  tields  as  the  shell  deforms.  Comparisons  between  the 
distorted  and  undistorted  elements  provides  dus  impression.  One  other  feature  is  that 
vibrational  characteristics  of  the  system  are  evident  by  examining  a  sequence  of  meshes 
separated  in  time.  Further,  gr{q)hicai  representations  provide  the  structural  arudyst  with 
information  on  which  criteria  can  be  based  or  a  method  to  focus-in  on  behaviors  that 
may  require  further  examination. 

7.2  Conclusions 

Numerically  stable  solutions  for  the  period  following  shell  collapse  are  always 
obtainable  with  the  iterative  methods  used  in  the  DSHELL  program.  The  pre-collapse 
time  step  selection  criteria  advocated  by  Greer  was  effective  for  all  cylindrical  shell 
configurations  considered.  The  post-coll{q)se  time  stq)  selection  criteria  considered  in  this 
research  was  effective  for  all  Donnell  shells,  but  not  for  all  SLR  shells.  The  SLR 
formulation  is  more  sensitive  to  time  step  choices  for  the  deeper  shells.  For  h/c  ratios  of 
0.273  (in  this  case,  the  6  in.  radius  shell)  and  greater,  successively  smaller  time  steps 
were  required  for  munerically  stable  results  during  the  coU^ed  shell  stages.  The  SLR 
formulation  carries  substantially  many  more  tenns  than  the  Donnell  formulation  and  is 
therefore  prone  to  a  combination  of  the  three  basic  numerical  errors,  namely  round-off, 
truncation,  and  propagated  errors  (16). 

For  predicting  the  shell  behavicn:  all  the  way  up  to  the  near  collapsed  configurations. 
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the  modified  Donnell  and  SLR  formulations  produce  essentially  the  same  results.  This 
is  an  important  result  because  the  Donnell  formulation  is  substantially  faster  than  the  SLR 
version.  Therefore,  the  Donnell  model  can  be  used  when  lai;ge  deflections  and  rotations 
are  not  a  prime  consideration.  However  in  predicting  the  behavitn*  at  collapse  and 
beyond,  the  two  models  produce  differing  results  for  the  deq)er  shell  dq>ths.  Shell  depth- 
to-chord  ratios,  /i/c,  less  than  or  i4)proximately  equal  to  0.276  (in  this  case,  the  8  in. 
radius  cylindrical  shell)  are  perhaps  the  limit  of  the  predictive  capability  of  the  SLR 
formulation;  this  is  still  a  substantial  number  since  it  implies  rotations  of  almost  60°  and 
a  displacement  that  is  50  times  the  thickness.  Although  it  is  recognized  that  plasticity  and 
material  failure  in  the  composite  would  occur  at  these  rotations  and  deflections,  it  does 
show  that  the  SLR  model  has  the  potential  to  accurately  predict  the  response  of  real 
systems,  if  plasticity  and  material  failure  were  incorporated  into  the  program,  since  the 
SLR  model  predicts  theoretical  behavior  that  cannot  be  achieved  by  current  composite 
materials.  At  these  hic  ratios,  the  Donnell  model  predicts  odd  shell  configurations  during 
the  transient  phase  and  much  larger  deflections  after  the  transient  phase  when  compared 
to  the  SLR  model.  The  Donnell  model  is  therefore  suspect  at  these  levels  of 
displacements  and  rotations.  It  is  not  clear  if  the  Donnell  model  is  suitable  for  predicting 
post  collapse  behavior  for  the  next  lowest  hIc  ratio  studied  (0.183)  since  the  shell  did  not 
collapse  as  expected.  This  phenomenon  shows  that  Donnell  models  tend  to  have  higher 
effective  flexural  stiffnesses,  as  one  might  expect  since  many  higher-order  rotatiorud  terms 
are  absent  from  the  formulation. 

For  shells  with  h/c  ratios  between  0. 183  and  0.076,  the  Donnell  and  SLR  formulations 
produce  similar  results  with  some  differences  during  the  transient  phase.  Therefore,  if  the 
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analyst  is  interested  in  long-term  behavior,  the  Donnell  model  is  sufficient  and  more 
efficient  from  a  computational  standpoint.  If  the  relatively  short  dynamic  phase  (which 
may  include  significant  backlash  for  the  structure)  is  an  important  consideration,  the  SLR 
model  is  recommended.  For  shells  with  hic  ratios  less  than  0.076,  including  plates,  the 
behaviors  predicted  by  the  two  models  are  similar.  Composite  plates  in  particular  had  no 
differences  in  the  center  node  displacement  curves.  Therefore  for  computational 
efficiency,  the  Donnell  model  should  always  be  used  when  very  shallow  shells  or  plates 
are  analyzed.  One  may  also  conclude  that  cylindrical  shells  with  hIc  ratios  less  than  or 
equal  to  0.076  is  one  defrnition  of  a  shallow  shell. 

Shells  that  do  not  colliq)se  have  the  potential  to  exhibit  chaotic  behavior.  The 
geometric  nonlinear  features  of  the  Donnell  and  SLR  formulations  makes  this  possible. 
Since  Greer  noted  harmonic  motion  with  loads  that  were  somewhat  less  than  the  collapse 
load  (about  70%  of  the  collapse  load),  loads  much  closer  to  the  coll^se  load  affect  the 
nonlinear  terms  to  the  extent  that  chaotic  motion  will  be  evident.  Transversely  loaded 
plates  do  not  collapse,  but  the  nonlinear  effects  will  result  in  the  chaotic  behavior  of 
composite  plates  as  well.  Moreover,  loading  rate  also  influences  the  nature  of  the  chaotic 
motion.  For  the  composite  plate,  high  loading  rates  results  in  erratic  phase-plane 
trajectories  whereas  gradual  loading  rates  facilitates  the  influence  of  the  strange  attractor. 

7.3  Recommendations 

The  graphical  displays  showed  dynamic  behaviors  not  noted  in  previous  research. 
These  behaviors  mainly  included  large  secondary  motions  during  the  transient  phase  and 
significant  in-plane  vibrational  motions.  These  observations  were  qualitative  in  nature. 
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Further  woik  can  be  done  to  characterize  diese  behaviors  in  a  quantitative  manner 
using  some  of  the  methods  used  in  this  work. 

In  order  to  more  fully  explore  the  issues  of  chaotic  behavior,  much  m(»e  data  is 
needed  than  was  generated  by  this  research.  Poincarg  maps  and  Lyapunov  estimations 
were  used  to  a  very  limited  extent  in  this  work  and  the  results  were  not  conclusive  due 
to  the  lack  of  data.  It  would  also  be  desirable  to  include  damping  effects  in  the  hope  of 
utilizing  other  methodologies  available  in  analyzing  chaotic  behavior  (damping  had  been 
used  in  earlier  research  with  DSHELL  (2)  but  to  a  very  limited  extent).  Computational 
efficiency  is  another  requirement.  Aside  from  improving  the  present  computer  hardware, 
one  suggestion  would  be  to  reduce  the  complexity  of  the  systems  to  be  analyzed.  Arches, 
for  example,  are  less  complex  than  cylindrical  shells  but  would  exhibit  many  of  the  same 
pre-  and  post-collapsed  characteristics;  computational  runs  should  produce  signiflcantly 
more  information  in  comparatively  shorter  times.  Another  possibility  is  simply  include 
damping  to  the  flat  plates  studied.  It  was  noted  that  the  Von  Karman  model  produced  the 
same  results  as  the  SLR  model.  The  Von  Karman  model  required  far  less  computational 
time  than  the  SLR  models.  With  more  data  generation,  more  definitive  results  with 
Ly^unov  exponents  and  Poincar6  m^  could  achieved.  Anotl^  quantitative  technique 
that  could  be  explored  is  the  calculation  of  fractal  dimensions.  Other  graphical  techniques 
such  as  logistic  maps,  circle  m^s,  horseshoe  maps,  etc.,  could  also  be  considered. 
Finally,  periodic  load  functions  should  be  considered  because  most  of  the  research  on 
chaotic  behavior  of  structural  systems  use  this  feature.  For  example.  Moon  (10: 133)  has 
shown  that  a  periodic  driving  force  facilities  the  identification  of  fractal-like  patterns. 
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APPENDIX  A  Effect  of  Radius  on  Flexural  Stiffness 


For  cylindrical  shells  with  the  same  cuivilinear  dimensions  (curvilinear  width,  length,  and 
thickness),  it  is  intuitive  that  a  shallower  shell  would  collapse  at  a  lower  transvere  loading 
than  a  deeper  shell.  Therefore,  the  effective  bending  stiffness  of  a  shell  is  affected  by  the 
radius.  To  illustrate  this  point,  consider  a  thin,  isotrc^ic  arch  with  a  constant  cross 
section: 


where  W  is  the  qjplied  load,  and  Afj,  H,  and  F,  are  the  moment  and  reaction  forces 
respectively.  In  order  to  simplify  the  problem  and  gain  a  general  idea  of  the  stress 
distributions,  the  arch  is  assumed  to  remain  rigid.  From  Roark  (18:163),  the  reaction 
forces  and  moment  is  given  by: 
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where  c  =  cos  0,  j  =  sin  0,  n  =  sin  (}>,  c  =  cos  (|»,  and  a  =  HAR.  I  and  A  are  the  moment 
of  inertia  and  aiea,  respectively,  of  the  cross  section.  /?  is  the  circular  radius  of  the  arch. 
For  this  example,  A  is  picked  as  a  square  cross  sectional  area  measuring  0.12"  x  0.12". 
The  applied  load  W  is  selected  as  100  lbs  and  is  ^plied  to  the  top  center  of  the  arch 
which  is  defined  as  the  zero  coordinate  position.  The  arch  length  is  also  held  constant 
at  12  in.  Using  the  above  equations  aiKl  resolving  H  and  V^  into  stresses  results  in 
Figures  A.2,  A.3  and  A.4.  These  figures  compare  the  stresses  and  moment  distributions 
as  a  function  of  the  arc  length.  Also  note  that  the  sign  convention  chosen  has  the 
compressive  stresses  as  positive  quantities. 

The  graphs  show  that  the  smaller  the  radius,  the  higher  the  normal  compressive  stresses. 
However,  the  stress  due  to  the  moments  must  be  taken  into  account.  The  stress  due  to 
moment  is; 
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Figure  A.4.  Radial  shear  stress  distribution  acting  along  an  arch 


where  Af,  is  moment  per  unit  length.  The  total  membrane  stress  is: 

A.5 

'  t  j’/12 

where  N,  is  normal  stress  per  unit  length.  Combining  the  results  of  Figures  A.2.  and  A.3. 
to  find  the  algebraic  sum  of  the  streeses  over  the  cross  sectional  area  results  in  Figure 
A.5.  which  shows  the  average  stress  distribution  along  the  arch.  The  figure  shows  that 
the  shallower  beam  has  generally  higher  levels  of  compressive  stresses.  Compressive 
membrane  forces  tend  to  reduce  the  bending  stifhiess  of  a  structure.  At  some  particular 
load  value,  the  compressive  stresses  reduce  the  bending  stiffiiess  to  zero  and  buckling 
occurs  [  :429].  For  more  complex  cases  such  as  a  thin  composite  cylindrical  shell,  the 
two  dimensionality  of  the  system  would  have  to  be  considered  as  well.  However,  this 
simple  example  implies  the  general  relationship  between  shell  radius  and  die  effective 
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bending  stiffness  for  thin  shells. 


Figure  Average  membrane  stress  distribution  along  a  thin  arch. 
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